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PREFACE 


Satellite-baaed altlmetrlc data taken by GEOS-3 and SEASAT over the Black Sea and Caspian 
Sea have been analyzed and a least-squares collocation technique used to predict the geoid 
undulations on a .25° x .25° grid and to transform these geoid undulations to free air gravity 
anomalies. The programs used to process and analyze these altimeter data are described and 
documented herein. A brief user’s guide is presented for each program that summarizes the 
purpose for, input to, and output from the code. This is followed by a complete listing of well 
commented code in the FORTRAN-77 language. These programs read, organize, and plot profiles 
and ground tracks of the altimeter data. They also identify crossovers between altimeter 
passes, determine and remove biases from the data, sort, and grid the corrected altimetric sea 
surface heights. Other programs generate contour plots of the grldded data, determine the 
gravity and geoid auto- and cross-covariance functions, and use these functions in a collocation 
estimate of gravity anomalies. 


lii 


PRECEDING PAGE BLANK NOT FILMED 



TABLE OF CONTENTS 


Preface m 

I INTRODUCTION : 

n USER'S GUIDES 2 

PRTDATA 2 

PLTGRP 4 

GRNTRK 5 

XOVERO 6 

XOVER 7 

PLTAEP 9 

SORT 10 

WGTAV.G ..!... 12 

CONTOUR 13 

GEOID 14 

EMPCQV 16 

GRAVEN ‘ ‘ ' 18 

JORDAN ‘ ‘ 21 

III REFERENCES 22 

IV LISTINGS 23 

PROGRAM PRTDATA 23 

PROGRAM PLTGRP 27 

PROGRAM GRNTRK 31 

PROGRAM XOVERO 35 

PROGRAM XOVER 44 

PROGRAM PLTAEP 60 

PROGRAM SORT 66 

PROGRAM WGTAVG | 71 

PROGRAM CONTOUR 76 

PROGRAM GEOID '''79 

PROGRAM EMPCOV 88 

PROGRAM GRAVEN. 91 

PROGRAM JORDAN. ’ i 07 


iv 



I INTRODUCTION 


The programs described in this document have been developed to process GEODYN-formatted 
satellite altimeter data, and to apply the processed results to predict geoid undulations and 
gravity anomalies of inland sea areas (Au et ai, 1988). These programs are written in standard 
FORTRAN 77 and are designed to run on the NSESCC IBM 3081 (MVS) computer. Because of 
the experimental nature of these programs, they are tailored to the geographical area analyzed. 
The attached program listings are customized for processing the altimeter data over the Black 
Sea. Users interested in the Caspian Sea data are expected to modify each program, although 
the required modifications are generally minor. Program control parameters are defined in the 
programs via PARAMETER statements and/or DATA statements. Other auxiliary parameters, 
such as labels, are ’hard-wired’ into the programs. Large data files are read in or written out 
through different input or output units. The program listings of these programs are 
accompanied by sample IBM Job Control Language ( JCL) Images. Familiarity with IBM JCL and 
the TEMPLATE graphic package is assumed. 
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n USER'S GUIDES 


PRTDATA 


This program reads a GEODYN-formatted altimeter data file, groups together the data records 
that belong to a satellite pass and extracts relevant geodetic and auxiliary information from 
each satellite pass. The process is repeated for all the passes in the input data file. For each 
satellite pass, the starting time in Julian day, the satellite identification number, the duration 
in seconds, the direction of the data arc (ascending or descending in latitude) and the number 
of data records associated with the pass are printed. For each data record, the corresponding 
record number in the file, the fraction of day in which the measurement Is taken, the absolute 
altimeter observation (meters), the ground track geodetic latitude and longitude and the observed 
surface elevation (meters) are printed. If a reference geoid undulation surface is supplied to 
the program, a comparison between the reference geoid undulation and the measured one is 
given for each data record. 

User-provided Parameters: 

MAXREC Maximum number of records In the input altimeter data file. 

XMIN Starting longitude of the input reference geoid undulation surface. 

XMAX Ending longitude of the input reference geoid undulation surface. 

YMIN Starting latitude of the input reference geoid undulation surface. 

YMAX Ending latitude of the input reference geoid undulation surface. 

DX Increment in longitude of the Input geoid reference undulation surface. 

DY Increment in latitude of the Input geoid reference undulation surface. 

Input: 

Reference geoid undulation surface (Unit 5): 
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The input file should be organized as follows: 

starting geodetic latitude YMIN (unformatted). 

geoid undulations of associated longitude from XMIN to XMAX 

incremented by DX (6(2X,F10.3)). 


YMIN + DY. 

geoid undulations of associated longitude from XMIN to XMAX incremented by DX. 


ending geodetic latitude YMAX. 

geoid undulations of associated longitude from XMIN to XMAX Incremented by DX. 

Altimeter data file (Unit 8): This is a binary file of altimeter data written in GEODYN 
format. 


Output: 

Organized printout (Unit 6). 

Reference geoid undulation corresponding to each data record (Unit 9): Each output record 
contains latitude, longitude and reference geoid undulation (5X,2F9.3,45X,F9.3). This 
output file may be used in subsequent programs for the purpose of comparison between 
measured geoid undulations and reference values. 
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PLTGRP 


This program reads a GEODYN-formatted altimeter data Ole, groups together the data records 
that belong to a satellite pass and plots the surface elevation profile of each pass as a function 
of time (fraction of day past midnight). The process is repeated for all the passes in the input 
data file. A surface elevation profile will not be plotted if there are less than four records 
In the satellite pass. The graphic software package, TEMPLATE, is used with this program. 

Input: 

Altimeter data Hie (Unit 5): This is a binary file of altimeter data written in GEODYN 
format. 


Output: 

Summary (Unit 6): Number of records read in input data file and the number of surface 
elevation profiles plotted. 

Plots of surface elevation profiles of each satellite pass as a function of time. 
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GRNTRK 


This program reads a GEODYN-formatted altimeter data file, groups together the data records 
that belong to a satellite pass and plots the ground track onto a geographical map. The 
process Is repeated for all the passes In the input data flle. A graphic software package. 
TEMPLATE, is used in this program. 


User-provided 

MAXREC 

NMAX 

VLONG1 

VLONG2 

VLAT1 

VLAT2 


Parameters: 

Maximum number of records in the input altimeter data flle. 
Maximum number of data records associated with one satellite pass. 
Starting east longitude of the output geographical map. 

Ending east longitude of the output geographical map. 

Starting latitude of the output geographical map. 

Ending latitude of the output geographical map. 


Input: 

Altimeter data flle (Unit 5): This Is a binary file of altimeter data written In GEODYN 
format. 


Output: 

Summary (Unit 6): Number of records read In the input flle. number of GEOS-3 and 
SEAS AT ground tracks plotted. 

Ground tracks plot. 
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XOVERO 


This program reads a GEODYN-formatted altimeter data file, approximates the ground track of 
each satellite pass as a second-degree equation and determines potential crossovers of these 
satellite passes. An overview of satellite crossovers is provided by this program before 
crossover adjustments are performed in a subsequent program, XOVER. Crossovers that occur 
at data gaps can be edited out 

User-provided Parameters: 

MPASS Maximum number of satellite passes in the Input altimeter data file. 

XMAX Maximum number of data records associated with a satellite pass. 

SIGE Discriminating criterion of a data spike. Observed surface elevation that deviates 

from adjacent one by more than SIGE will be treated as noise and, thereby, 
will not be included in subsequent data analyses. 

Input: 

Altimeter data file (Unit 3): This is a binary file of altimeter data written in GEODYN 
format. 


Output: 

Summary (Unit 6): Possible satellite crossovers, including sateilite types, crossover 
locations. 
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XOVER 


Crossover adjustment is performed in this program. This program reads a GEODYN-formatted 
altimeter data file, approximates each satellite pass ground track by a second-degree equation 
and determines possible crossovers of these satellite passes. The geoid height at each 
crossover location for each pass is obtained by interpolating screened data in that pass. A 
weighted least-squares analysis is performed to determine biases for each pass so that the 
crossover restduais are minimized, holding one pass fixed as the reference pass. All satellite 
passes are then rectified to a common reference geoid through the chosen reference pass. The 
output includes the RMS of crossover residual before and after orbit bias adjustment, a 
GEODYN-formatted file with the original observed surface elevation replaced by the adjusted 
one, and an error covariance matrix of the bias adjustment. 

User-provided Parameters: 

SIGMAG Assumed noise for GEOS-3 orbit. 

SIGMAS Assumed noise for SEASAT orbit. 

XDMAX Maximum crossover residual for that crossover to be considered. 

ETOP Upper limit of observed surface elevation to be accepted. 

EBOT Lower limit of observed surface elevation to be considered. 

MAXOVR Maximum number of crossovers. 

MPASS Maximum number of satellite passes in the input altimeter data file. 

NMAX Maximum number of data records associated with a satellite pass. 

SIGE Discrimination criterion of a data spike. Observed surface eievation that deviates 

from adjacent one by more than SIGE will be treated as noise and, thereby, 
will not be included in subsequent data analyses. 

Input: 

Reference geoid undulation model (Unit 2): Each record in the file contains latitude, 
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longitude and reference geoid undulation (5X.2F9.3.45X.F9.3). This is the output file 
from program PRTDATA. 

Altimeter data file (Unit 3): This is a binary file of altimeter data written in GEODYN 
format 

Crossovers to be edited out (Unit 4): This file contains records of two passes that are to 
be excluded from crossover adjustment Each record contains three variables, JC, IP, 
and JP, where IP and JP are the pass number of the two passes that form the 
crossover, and JC is an editing criterion. If JC is set to zero, that crossover will be 
excluded from crossover adjustment The format of the record is (11,2110). 


Output: 

Adjusted altimeter data file (Unit 8): This is a binary file written In GEODYN format with 
the observed surface elevation replaced by the adjusted one. 

Error covariance matrix of the bias adjustment (Unit 9): The format of this output matrix 
is (4D20. 13). 
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PLTAEP 


This program reads an adjusted altimeter data Ole in GEODYN format and a file containing 
reference geoid undulation that is in one-to-one correspondence with the adjusted altimeter 
data, and plots the adjusted elevation profile of each satellite pass. The reference geoid is 
also plotted onto each elevation profile for comparison. A graphic software package, 
TEMPLATE, is used in this program. 

User-provided Parameters: 

MAXPLT Maximum number of plots to be performed. 

ETOP Upper limit of observed surface elevation to be accepted. 

EBOT Lower limit of observed surface elevation to be considered. 

NMAX Maximum number of data records associated with a satellite pass. 

NINT Number of intervals used in a smoothing routine to plot out a smoothed elevation 

profile of a satellite pass. 

SIGE Discriminating criterion of a data Spike. Observed surface elevation that deviates 

from adjacent one by more than SIGE will be treated as noise and, thereby, 
will not be included in subsequent data analyses. 

SIGT Not used. 

Input: 

Reference geoid undulation model (Unit 3): Each record in the file contains latitude, 

longitude and reference geoid undulation (5X,2F9.3,45X,F9.3). This is the output file 
from program PRTDATA. 

Altimeter data file (Unit 4): This is a binary file of altimeter data written in GEODYN 
format. 

Output: 

Plots of altimeter surface elevation profiles adjusted for orbit biases. A reference geoid is 
superimposed on each profile for comparison. 
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SORT 


This program reads an adjusted altimeter data file and determines with which grid points a 
data record should be associated, according to a specified cap size. Each GEODYN data record 
is checked for grid point association. For each grid point associated with a record, an output 
record is written with latitude, longitude, adjusted surface eievation, grid points identification 
number, grid point counter, associated pass number, and standard deviation of that data point. 
A sort/merge process (TSORT In MVS system) can be performed to group together ail the 
altimeter data points that are associated to each grid point. 

User-provided Parameters: 

ETOP Upper limit of adjusted surface elevation to be accepted. 

EBOT Lower limit of adjusted surface eievation to be considered. 

SIGE Discriminating criterion of a data spike. Observed surface eievation that deviates 

from adjacent one by more than SIGE will be treated as noise and, thereby, 
will not be included in subsequent data analyses. 

MIREC Maximum number of data records associated with a satellite pass. 

XMIN Starting longitude of the chosen grid. 

XMAX Ending longitude of the chosen grid. 

YMIN Starting latitude of the chosen grid. 

YMAX Ending latitude of the chosen grid. 

DX Increment in longitude of the chosen grid. 

DY Increment in latitude of the chosen grid. 

CAP Cap size within which altimeter data is considered. 

Input: 

Altimeter data file (Unit 5): This is a binary file of adjusted altimeter data written in 
GEODYN format. 
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Output: 


Summary (Unit 6): Data partitioning information and data distribution map. 

Altimeter data partition (lie (Unit 8): Each record in this binary file contains the 

longitude, latitude, adjusted surface elevation, grid identification index for 
longitude, grid identification index for latitude, satellite pass number, and 
standard deviation of the adjusted surface elevation. In order to be used in 
subsequent Program GEOID, this output file should be sorted in ascending order 
of the grid longitude index and grid latitude index. 
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WGTAVG 


This program reads an adjusted altimeter data file written in GEODYN format, sorts each data 
record into a .25° x .25° grid, and determines a weighted average of the data points at each 
grid point. 

User-provided Parameters: 

ETOP Upper limit of adjusted surface elevation to be accepted. 

EBOT Lower limit of adjusted surface elevation to be considered. 

SIGE Discrimination criterion of a data spike. Observed surface elevation that deviates 

from adjacent one by more than SIGE will be treated as noise and, thereby, 
will not be included in subsequent data analyses. 

MIREC Maximum number of data records associated with a satellite pass. 

XMIN Starting longitude of the chosen grid. 

XMAX Ending longitude of the chosen grid. 

YMIN Starting latitude of the chosen grid. 

YMAX Ending latitude of the chosen grid. 

DX Increment in longitude of the chosen grid. 

DY Increment in latitude of the chosen grid. 

MXP Maximum number of data points associated with a grid point. 

Input: 

Altimeter data file (Unit 5): This is a binary file of adjusted altimeter data written in 
GEODYN format. 

Output: 

Summary (Unit 6): Data partitioning information, weighted averages and data distribution 
map. 

Gridded weighted -averages (Unit 8): The output file contains sequences of grid longitude, 
grid latitude and surface elevation. The format of the output file is (3(2F7.2,F12.4). 
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CONTOUR 


This program reads a grldded data file and plots the corresponding contour. A geographical 
map defined within the chosen grid boundary is superimposed on the contour. Any grid point 
that has no physical value should be masked with the value of -1000.0 and It will not be 
contoured. A graphic software package, TEMPLATE, is used In this program. 


User-provided 

BIAS 

XMIN 

XMAX 

YMIN 

YMAX 

DX 

DY 

CHGH 

CLOW 

CINC 

CLAB 


Parameters: 

Overall bias that Is to be subtracted of T all the data points. 
Starting longitude of the chosen grid. 

Ending longitude of the chosen grid. 

Starting latitude of the chosen grid. 

Ending latitude of the chosen grid. 

Increment in longitude of the chosen grid. 

Increment In latitude of the chosen grid. 

Upper limit of contour. 

Lower limit of contour. 

Increment interval of contour. 

Frequency of labelling the contour lines. 


Input: 

Grldded data file (Unit 5): The output file contains sequences of grid longitude, grid 
latitude and surface elevation. The format of the output file is (3(2F7,2,F 12.4). 


Output: 

A contour map of the gridded data. 
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GEOID 


This program predicts and grids geoid undulations using a collocation technique. The input 
dataset is a sorted version, according to grid longitude and latitude indices, of the output file 
from Program SORT. 


User-provided Parameters: 

MPASS Number of satellite passes considered in the adjusted altimeter dataset. 

NREFP Relative position number of the reference satellite pass in the aitimeter dataset. 

KUTOFF Minimum number of data points associated with a grid point to be considered 

for geoid prediction. 

LCTV Set to 1. 

MTR Maximum number of data records associated with a grid point. 

XMIN Starting longitude of the chosen grid. 

XMAX Ending longitude of the chosen grid. 

YMIN Starting latitude of the chosen grid. 

YMAX Ending latitude of the chosen grid. 

DX Increment in longitude of the chosen grid. 

DY Increment in latitude of the chosen grid. 

LCT Number of intervals in the given covariance function. 

Input: 

Covariance table (Unit 5): Each record in the covariance table contains distance in 
degree, geoid -to-geoid value, geoid-to-gravity value, and gravity-to-gravity value. 

The format of the input file is ( 1X,D 10.4.3D20. 13). 

Sorted altimeter data file (Unit 8): This is a sorted version of the binary output file 

from Program SORT. The sorting should be performed in ascending order of grid 
longitude index and then of grid latitude index. 

Error covariance matrix of orbit adjustment (Unit 11): This is the output error 
covariance matrix from Program XOVER. 
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Output: 


Gridded geoid undulations (Unit 9): The output file contains sequences of grid 

longitude, grid latitude and surface elevation. The format of the output file is 
(3(2F7.2,F 12.4). 

Gridded variance of predictions (Unit 10): The output file contains sequences of grid 
longitude, grid latitude and variance in cm. The format of the output file is 
(3(2F7.2,F12.4). 
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EMPCOV 


This program determines a set of residual local empirical covariance functions based on a 
convolution technique. Reference models as well as gridded data of both geoid undulations and 
gravity anomalies must be supplied to the program 

User-provided Parameters: 

XMIN Starting longitude of the chosen grid. 

XMAX Ending longitude of the chosen grid. 

YMIN Starting latitude of the chosen grid. 

YMAX Ending latitude of the chosen grid. 

DXY Increment in longitude and latitude of the chosen grid. 

Input: 

Reference geoid undulation surface (Unit 7): 

The input file should be organized as follows: 
starting geodetic latitude YMIN (unformatted). 

geoid undulations of associated longitude from XMIN to XMAX incremented 
by DXY (6(2X.F 10.3)). 

YMIN + DXY. 

geoid undulations of associated longitude from XMIN to XMAX Incremented by DX. 
ending geodetic latitude YMAX. 

geoid undulations of associated longitude from XMIN to XMAX incremented by DX. 
Reference gravity anomalies surface (Unit 9): 

Input format similar to Unit 7. 

Gridded geoid undulations (Unit 8): The input file contains sequences of grid longitude, 
grid latitude and surface elevation. The format of the input (lie is (3(2F7.2,F12.4). 

Gridded gravity anomalies (Unit 10): The input (lie contains sequences of grid 
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longitude, grid latitude and gravity anomalies. The format of the input file Is 
(3(2F7.2,F12.4). 


Output: 

Covariance table (Unit 11): Each record in the covariance table contains distance in 
degree, geoid-to-geoid value, geoid-to-gravity value, and gravlty-to-gravlty value. 
The format of the output file is ( 1X,D 10.4.3D20. 13). 
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GRAVEN 


This program transforms gridded geoid undulation data into gravity anomalies. Three reference 
models. Rapp‘s 36 x 36. 180 x 180. and 300 x 300 gravity models (Rapp. 1986) can be used as 
reference surfaces in this program. An iterative process to refine the reference surfaces can 
also be performed by this program. 


User-provided Parameters: 


MODEL 

ANOISE 

ITM 

CAP 

IPRT 

CNVRGE 

SCX 

KUTOFF 

MTR 

XMIN 

XMAX 

YMIN 

YMAX 

XG1 

XG2 

YG1 

YG2 


Reference model to be used. It can be 36. 180 or 300. which corresponds to 
Rapp’s 36 x 36. 180 x 180. and 300 x 300 models. 

Assumed noise of the geoid undulation data. 

Maximum number of iterations to be carried out. When there is no iteration to 
be performed. ITM is set to zero. 

Integration cap size used. 

Option for printouts. Set IPRT to zero for minimum output. 

Convergence criterion for iterative process. 

Set to 1. 

Minimum number of data points associated with a grid point to be considered for 
gravity transformation. 

Maximum number of data records associated with a grid point. 

Starting longitude of the chosen grid. 

Ending longitude of the chosen grid. 

Starting latitude of the chosen grid. 

Ending latitude of the chosen grid. 

Set equal to XMIN. 

Set equal to XMAX. 

Set equal to YMIN. 

Set equal to YMAX. 
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DX 


Increment in longitude of the chosen grid. 
Increment in latitude of the chosen grid. 


DY 


Input: 

Covariance table (Unit 5): Each record in the covariance table contains distance in 
degree, geoid-to-geoid value, geoid-to-gravity value, and gravity-to-gravlty value. 

The format of the input (lie is ( IX, D L0.4,3D20. 13). 

Gridded geoid undulations (Unit 8): The input file contains sequences of grid longitude, 
grid latitude and surface elevation. The format of the input file is (3(2F7.2,F12.4). 

Gridded variance of geoid prediction (Unit 9): The Input file contains sequences of grid 
longitude, grid latitude and variance in cm. The format of the input file is 
(3(2F7.2*F 12.4). 

Reference geoid undulation surface. Rapp’s 180 x 180 model (Unit 10): 

The input file should be organized as follows: 
starting geodetic latitude YMIN (unformatted). 

geoid undulations of associated longitude from XMIN to XMAX incremented 
by DX (6(2X,F10.3)). 


YMIN + DY. 

geoid undulations of associated longitude from XMIN to XMAX incremented by DX. 
ending geodetic latitude YMAX. 

geoid undulations of associated longitude from XMIN to XMAX incremented by DX. 


Reference gravity anomalies surface, Rapp's 180 x 180 model (Unit 11): 
Input format similar to Unit 10. 


Reference geoid undulation surface, Rapp’s 36 x 36 model (Unit 12): 
Input format similar to Unit 10. 

Reference gravity anomalies surface, Rapp s 36 x 36 model (Unit 13): 
Input format similar to Unit 10. 

Reference geoid undulation surface, Rapp’s 300 x 300 model (Unit 14): 
Input format similar to Unit 10. 
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Reference gravity anomalies surface, Rapp's 300 x 300 model (Unit 15): 


Input format similar to Unit 10. 


Mask of inland sea surface (Unit 16): This is an M x N mask of the inland sea surface, 
where M and N are the numbers of longitude and latitude grid points. A unitary 
value ( 1) indicates that a prediction and transformation are to be performed at 
that grid point. A zero (0) indicates that no prediction or transformation is to 
take place at that grid point. 


Output: 

Gridded gravity anomalies (Unit 17): The output file contains sequences of grid 

longitude, grid latitude and gravity anomalies. The format of the output flic i9 
(3(2F7.2,F12.4). 

Gridded variance of gravity transformation (Unit 18): The output file contains 

sequences of grid longitude, grid latitude and variance in mgal. The format of the 
output file is (3(2F7.2,F12.4). 

Gridded geoid undulations (Unit 19): The output file contains sequences of grid 

longitude, grid latitude and geoid undulations. The format of the output Hie is 
(3(2F7.2,F12.4). 

Gridded variance of geoid prediction (Unit 20): The output file contains sequences of 
grid longitude, grid latitude and variance In cm. The format of the output file is 
(3(2F7.2,F 12.4). 

Gridded difference between gravity anomalies at the end of the Iterative process and 
the reference surface used (Unit 21): The output file contains sequences of grid 
longitude, grid latitude and gravity anomalies. The format of the output file is 
(3(2F7.2,F12.4). 

Covariance table (Unit 22): This is a set of residual local empirical covariance 
functions obtained at the end of the iterative process. Each record in the 
covariance table contains distance In degree, geoid-to-geold value, geoid-to-gravity 
value, and gravlty-to-gravity value. The format of the output file is 
( IX, D 10.4, 3D20. 13). 
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JORDAN 


This program determines a set of theoretical covariance functions according to Jordan’s 
formation (Jordan, 1972). 


User-provided Parameters: 

SN Variance of geoid undulations at prediction point. 

SG Variance of gravity anomalies at prediction point. 

D Correlation distance. 

DR Increment in distance in degrees. 

RMIN Starting value of the distance In degrees. 

RMAX Ending value of the distance In degrees. 


Output: 

Covariance table (Unit 8): Each record in the covariance table contains distance in 
degree, geoid-to-geoid value, geoid-to-gravity value, and gravity- to-gravity value. 
The format of the output file is ( 1X,D 10,4, 3D20. 13). 
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IV LISTINGS 


A listing of program prtdata 


/ /ZMAYABSR JOB ( G0109 , 360 , 2 ) , AYAU, TIME- ( 0 , 30 ) , CLASS-O, MSGCLASS-X 
/♦JOBPARM LINES-60 
// EXEC FORTVC 
//SYSIN DO * 

C 

C FORMAT OF GEOS-3/SEASAT ALTIMETER DATA 


c 

VARIABLE 

TYPE 

c 



c 

★ 

14 

c 


12 

c 


12 

c 


14 

c 


14 

c 

* 

14 

c 

* 

R8 

c 

★ 

R8 

c 

* 

14 

c 

♦ 

14 

c 

* 

R4 

c 


14 

c 


14 

c 


14 

c 


R4 

c 


14 

c 

* 

14 

c 

★ 

14 

c 


14 

c 


14 

c 


12 

c 


12 

c 


12 

c 


12 

c 


12 

c 


12 

c 


12 

c 


12 

c 


12 

c 


12 

c 


12 

c 


12 


C 

C 


c 

c 


DESCRIPTION 

SATELITE ID 
MEASUREMENT TYPE < 

TIME SYSTEM ( NM ) 

STATION NUMBER 

PREPROCESSING INDICATORS 

MODIFIED JULIAN DATE OF OBSERVATION 


42* OVER LAND , 43 » OVER WATER 


(GMT) 

( METERS ) 

( IE-6 DEGREES) 
( IE-6 DEGREES ) 
( METERS ) 

(MM) 

(GEODYN VOL 3) 
(MM) 


FRACTION OF DAY PAST MIDNIGHT 
ALTIMETER OBSERVATION 
SATELLITE GEODETIC LATITUDE 
SATELLITE EAST LONGITUDE 
MEASUREMENT STANDARD DEVIATION 
NET INSTRUMENT CORRECTION 
METEOROLOGICAL DATA WORD 
NET MEDIA CORRECTIONS 

GEOID HEIGHT ABOVE REFERENCE ELLIPSOID (METERS) 

NET OCEAN DYNAMIC CORRECTIONS (MM) 

INDICATED SURFACE ELEVATION (MM) 

S/C REVOLUTION NUMBER 

MEAN SEA SURFACE ELEVATION (MARSH/MARTIN '81 (MM)) 
DOD REFERENCE RADIAL ORBIT DIFFERENCE (MM) 

H 1/3 (CM) 

AGC ( DB ) 

WIND SPEED (CM/SEC) 

SURFACE ELEVATION PREPROCESSING WORD 
DRY TROPOSPHERIC CORRECTION (MM) 

FNOC WET TROPOSPHERIC CORRECTION (MM) 

SMMR WET TROPOSPHERIC CORRECTION (MM) 

IONOSPHERIC CORRECTION (MM) 

BAROTROPIC DYNAMIC SEA SURFACE CORRECTION (MM) 

(MM) 

(MM) 

(MM) 


SOLID EARTH TIDE 
SCHWIDERSKI OCEAN TIDE 
PARKE OCEAN TIDE 


PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 


INTEGER*2 
INTEGER* 4 
REAL *8 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
REAL *4 


MAXREC 
XMIN - 
XMAX * 
YMIN = 
YMAX » 
DX = 0 
DY = 0 
NLON => 
NLAT = 
IC 
IWK 


10000 

) 


) 


- 26.5 

* 42.5 ) 

=* 40.0 ) 

* 48.0 ) 

0.25 ) 

0.25 ) 

* ( (XMAX-XMIN) /DX + 1 ) ) 

* ( ( YMAX- YMIN) /DY + 1 ) ) 
NLON ) 

2* NLON* NLAT + 2 * ( NLON+NLAT ) 


) 


12,13,117,118,119,120,121,122,123,124.125,126,127,128 
II, 14, 15, 16, 17, 18, 19, 110, 111, 112, 113, 114, 115, 116 
R1 , R2 
R3,R4 

LON (NLON), LAT(NLAT), GRID( NLON, NLAT) 

C( 2, NLON, 2, NLAT) , WK(IWK) 

GLAT, ELON, HSS, GUND 

DUR, SLAT1 , SLAT2 , SLON1, SLON2 
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V ouu uuuo 


c 


REAL* 8 

INTEGER 

INTEGER 

INTEGER 

CHARACTER 

CHARACTER 

CHARACTER 

CHARACTER 

CHARACTER 


STIM1 , STIM2 

ING, INP, OUP, OUR 

IREC, LPP, NLINE, NPAGE , NREC 

IIS, I6S, IDUMl 

TITLE*20 /' BLACK SEA DATA '/ 
NS* 12 /' NORTH-SOUTH'/ 

SN*12 /' SOUTH-NORTH'/ 

EW* 12 /'SINGLE-POINT'/ 

DIR* 12 


D*™ / 5 /, INP / 8 /, OUP / 6 /, OUR / 9 / 

DATA K75 / 0 /, K78 / 0 / 

DATA LPP / 54 /, NLINE / 1 /, NPAGE / 1 /, NREC / 0 / 


GENERATE A GRID NET AND READ IN RAPP'S OSU300 180X180 GEOID 


LON ( 1 ) - XMIN 
DO 100 1-2, NLON 

LON ( I ) - LON(I-l) 
100 CONTINUE 

LAT< 1) - YMIN 
DO 110 J » 2, NLAT 
LAT(J) - LAT(J-l) 
110 CONTINUE 


+ DX 


+ DY 


DO 120 
READ < 
READ ( 
120 CONTINUE 


J - NLAT, 1, -1 
ING, * ) RL 

ING, 501 ) (GRID(I,J) ,1-1, NLON) 


DETERMINE A BICUBIC SPLINE TWO-DIMENSIONAL COEFFICIENT MATRIX 
CALL IBCCCU ( GRID, LON, NLON, LAT , NLAT , C, IC, WK, IER ) 


C 

C 


C 

C 

c 


START PROCESSING THE GEODYN DATA SET 


READ ( INP ) 


& 

6 


NREV - 114 
IIS - II 
I6S - 16 


? 4 ' 15 ' i6 »R1»R2,I7, 18, R3, 

19,110,111, R4 ,112,113,114,115,116,117.118 119 
120,121,122,123,124,125,126,127,128 '' 


SLAT1 
SLON1 
STIM1 - R1 
IF ( II. EQ 
IF ( II 
WRITE ( 
WRITE ( 


1.0E-6 
1 . 0E-6 


FLOAT { 17) 
FLOAT (18) 


REWIND INP 


7502701 ) K75 
EQ. 7806401 ) K78 
OUP, 601 ) TITLE, 
OUP, 602 ) 


- K75 + 1 
= K78 + 1 

16, NREV, II, NPAGE 


DO 111 MJ - 1, MAXREC 


READ ( INP, END 

& 

fc 

NREC - NREC + 1 
GLAT = 1 . OE-6 * 
ELON - 1 . OE-6 * 
HSS - 1.0E-3 * 


tq 2 ??h M 1 '5 2 '^' I4 ' I5 ' I6 ' R1 ' R2 ,I7.I8,R3, 

,1 ,1 , R ,ii2,ii3,n 4 .n 5 ,ii 6 ,n 7 ,ii8,xi 9 , 

120,121,122,123,124, 125,126,127,128 


FLOAT ( 17) 
FLOAT ( 18) 
FLOAT ( 113) 


CALL IBCEVL ( LON, NLON, LAT, NLAT, C,IC, ELON, GLAT, GUN D, IER ) 

WRITE < OUR, 660 ) GLAT, ELON, GUND 

IF ( 114. EQ. NREV ) THEN 
IF ( NLINE. LE. LPP ) THEN 

WRITE ( OUP, 603 ) NREC, R1,R2, GLAT, ELON, HSS, R3, GUND 
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NLINE - NLINE + 1 
IREC - IREC + 1 
ELSE 

NLINE « 1 

IREC - IREC + 1 

NPAGE - NPAGE + 1 

IF ( (I6-I6S). EQ. 1 ) R1 ■ R1 + 1.0 

WRITE ( OUP, 601 ) TITLE, 16, NREV, II, NPAGE 

WRITE < OUP, 602 ) 

WRITE ( OUP, 603 ) NREC , R1 , R2 , GLAT , ELON , HSS , R3 , GUND 
END IF 

SLAT2 - GLAT 
SLON2 - ELON 
STIM2 - R1 
ELSE 

IF ( SLAT2. EQ. 0.0. AND. SLON2 . EQ. 0.0 ) THEN 
DIR - EW 
DUR - 0.0 
ELSE 

IF ( (SLAT2-SLAT1) . LE. 0.0 ) THEN 
DIR - NS 
ELSE 

DIR - SN 
END IF 

DUR - 86400.0 * ( STIM2 - STIMl ) 

END IF 

WRITE ( OUP, 610 ) DIR, DUR, IREC 

IF ( II. EQ. 7502701 ) K75 - K75 + 1 

IF ( II. EQ. 7806401 ) K78 - K78 + 1 

NREV - 114 

IIS - II 

I6S - 16 

NLINE - 1 

IREC * 1 

NPAGE - NPAGE + 1 

WRITE ( OUP, 601 ) TITLE, 16, NREV, II, NPAGE 
WRITE ( OUP, 602 ) 

WRITE ( OUP, 603 ) NREC , R1 , R2 , GLAT, ELON, HSS , R3 , GUND 

SLAT1 - GLAT 

SLON1 - ELON 

STIMl - R1 

SLAT 2 - 0. 

SLON2 * 0. 

STIM2 » 0. 

END IF 

111 CONTINUE 

2000 CONTINUE 

IF ( SLAT2 . EQ. 0.0. AND. SLON2 . EQ. 0.0 ) THEN 
DIR » EW 
DUR =» 0.0 
ELSE 

IF ( ( SLAT2-SLAT1 ) . LE. 0.0 ) THEN 

DIR = NS 
ELSE 

DIR =■ SN 
END IF 

DUR = 86400.0 * ( STIM2 - STIMl ) 

END IF 

WRITE ( OUP, 610 ) DIR, DUR, IREC 
WRITE ( OUP, 611 ) K75, K78 

501 FORMAT ( (6(2X,F10.3) ) ) 

601 FORMAT ( ' 1 ’ , T2 , A2 0 , 8X , ' MOD . JULIAN DAY = ’ , 16 , 8X, ' PASS = ', 

«. 8X, 'SATELLITE ID = ' , 1 10 , T120 , ' PAGE ',14,/) 

602 FORMAT ( ' ',' RECORD # ',2X,' FRACTION OF DAY',3X, 

* ' ALTIMETER OBS . (M)',5X, ' S/C GEOD LAT ', 


4 

4 

& 

4 


' S/C EAST LONG ' 
' + ' / ' ' ,2X, ' 


ELEVATION (M)',5X, ' RAPP- 180 '/ 
- ,3X ' 


, ■' l uv t 

S0 V“ * MXIr 

604 FORMAT ( 15,110, 2F8 . 3, 4F10 . 6 ) 

610 FORMAT(lH ,/,T5,'THE APPROX. DURATION OP THIS ',A12,' PASS IS' 

,,, 4 _ F8. 2 , SEC. WITH ',16,' RECORDS.') IS ' 

611 roRMAT ( !H1,T20, ' # OF PASSES FOR SATELLITE GEOS-3 ISs' 16 // 

-, 4 _ _ T20, ' # OF PASSES FOR SATELLITE SEASAT IS*' IS I 

660 FORMAT< 5X, 2F9 . 3, 45X, F9 . 3 ) seasat IS. ,16) 


STOP 

END 

II* 

// EXEC LINKGOV, REGION. GO-5000K 
//SYSLIB DD DSN-SYS2.IMSLS,DISP-SHR 

DD dsn “ z MAYA.ALTIM.DATA(BRH180) ,DISP-SHR 
/ /GO. FT06F001 DD SYSOUT-* " K 

//GO.FTOSrOOl DD dsn-zmaya. black. data, disp-shr 
//* ***. **** ...* **** ***. **,* **** # 

If* BLACK SEA DATA 

II* 

H* **** .... ...» .... ..** „„ .... 

//* ** * 

//GO.FT09F001 DD DSN-ZMAYA. ALTIM. DATA ( BLKRAPP ), DISP-SHR 

// EXEC NOTirYTS 
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nnnnooooonnnoooooononooooonnononooon 


A listing Of PROGRAM PLTGRP 


/ /ZMAYABEP JOB ( GO 1 09 , 360 , 75 ) , AYAU, TIME»( 12 , 00 ) , CLASS-O, MSGCLASS-X 
/* JOBPARM LINES-60 
//TEMPLATE EXEC G38PLOT 
//rORT.SYSIN DD * 


FORMAT OF SEASAT ALTIMETER DATA 
VARIABLE TYPE DESCRIPTION 

SATELITE ID 
MEASUREMENT TYPE ( 42 
TIME SYSTEM ( NM ) 

STATION NUMBER 
PREPROCESSING INDICATORS 
MODIFIED JULIAN DATE OF OBSERVATION 
FRACTION OF DAY PAST MIDNIGHT (GMT) 
ALTIMETER OBSERVATION 
SATELLITE GEODETIC LATITUDE 
SATELLITE EAST LONGITUDE 


TYPE 

14 

12 

12 

14 

14 

14 

R8 

R8 

14 

14 

R4 

14 

14 

14 

R4 

14 

14 

14 

14 

14 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 


OVER LAND , 43 = OVER WATER ) 


( METERS ) 

( IE-6 DEGREES ) 
( IE-6 DEGREES) 


MEASUREMENT STANDARD DEVIATION (METERS) 


(MM) 

(GEODYN VOL 3) 
(MM) 


NET INSTRUMENT CORRECTION 
METEOROLOGICAL DATA WORD 
NET MEDIA CORRECTIONS 

GEOID HEIGHT ABOVE REFERENCE ELLIPSOID (METERS) 

NET OCEAN DYNAMIC CORRECTIONS (MM) 

INDICATED SURFACE ELEVATION (MM) 

S/C REVOLUTION NUMBER 

MEAN SEA SURFACE ELEVATION (MARSH/MARTIN '81 (MM)) 
DOD REFERENCE RADIAL ORBIT DIFFERENCE (MM) 

H 1/3 (CM) 

AGC ( DB ) 

WIND SPEED (CM/SEC) 

SURFACE ELEVATION PREPROCESSING WORD 
DRY TROPOSPHERIC CORRECTION (MM) 

FNOC WET TROPOSPHERIC CORRECTION (MM) 

SMMR WET TROPOSPHERIC CORRECTION (MM) 

IONOSPHERIC CORRECTION (MM) 

BAROTROPIC DYNAMIC SEA SURFACE CORRECTION (MM) 
SOLID EARTH TIDE (MM) 

SCHWIDERSKI OCEAN TIDE (MM) 

PARKE OCEAN TIDE (MM) 


C 


C 


INTEGER* 2 
INTEGER* 4 
REAL* 8 
REAL* 4 
REAL *4 
REAL *4 
INTEGER 
INTEGER 
CHARACTER 


12, 13, I 17, I 18, I 19, 120, 121, 122, 123, 124, 125, 126, 127, 128 
II, 14, 15,16,17, 18, 19, 110, 111,112,113,114, 115, 116 
R1,R2 
R3,R4 

GLAT, ELON, HSS 
X( 5000) ,Y( 5000) 

INP, IREC, LPP , NLINE, NPAGE, IPLOT, NREC, NREV, OUP 
IWRT , MAXPLT, LIST(6) 

TITLE*20 /' BLACK SEA DATA '/ 


DATA INP / 5 /, OUP / 6 / 

DATA LIST / 8, 110, 75, 0, 7, 10 / 

DATA LPP / 54 /, NLINE III, NPAGE III 
DATA IREC / 0 /, NREC 101, IPLOT 101 
DATA IWRT 101, MAXPLT / 99999 / 


READ ( INP ) I1,I2,I3,I4,I5,I6,R1,R2,I7,I8,R3, 

4 19, 110, 111, R4, 112, 113,114, 115,116, 117, 118, 119, 

4 120,121,122,123,124,125,126,127,128 

NREV =114 

IF ( IWRT. GT. 0 ) THEN 

WRITE ( OUP, 601 ) TITLE, NREV, NPAGE 
WRITE ( OUP, 602 ) 

END IF 

REWIND INP 

CALL USTART 

CALL UINQES <1.0, SUPRT ) 
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non 


IF ( SUPRT . NE .0.0) CALL UESCAP( 1 . 0,LIST, 6 . 0, IDUM,DUM) 


1000 READ 
★ 

* 

NREC 

GLAT 

ELON 

HSS 


( INP, END 


NREC + 
1. OE-6 
1.0E-6 
1.0E-3 


ra 2 r?2 > I1,I2,I3,I4,I5,I6,R1,R2,I7,I8,R3, 

FLOAT ( 17) 

FLOAT ( 18) 

FLOAT ( 113) 


IF ( II 4, EQ. NREV ) THEN 
IF ( NLINE. LE . LPP ) THEN 

IF (IWRT.GT.0) WRITE ( OUP , 603) 
NLINE - NLINE +1 
ELSE 

NLINE - 1 
NPAGE « NPAGE + 1 
IF ( IWRT. GT. 0 ) 

WRITE ( OUP, 601 


16 , R1 , R2 , GLAT, ELON , HSS , NREC 


WRITE { 

WRITE ( 

END IF 
END IF 
IREC - IREC 
RMJD - FLOAT (16) 
RNREV * FLOAT (114) 
X(IREC) - R1 
Y(IREC) * HSS 
ELSE 


OUP, 

OUP, 


+ 1 


602 

603 


THEN 

TITLE, NREV, NPAGE 


) 16 , RI , R2 , GLAT , ELON , HSS , NREC 


PLOT THE PREVIOUS ELEVATION PROFILE 


IF ( I PLOT . GT. MAXPLT ) GO TO 3000 

IF ( IREC. GT. 3. AND. IPLOT. LE . MAXPLT ) THEN 

I PLOT * IPLOT + 1 

CALL PLOT(X, Y, IREC, RMJD, RNREV) 

END IF 


NREV - 114 
NLINE - 1 
NPAGE - NPAGE + 1 
IREC - 1 


X< IREC) 

m 

Rl 


Y( IREC) 

m 

HSS 


IF ( IWRT 

. GT. 

o ) 
601 

WRITE 

( 

OUP, 

WRITE 

{ 

OUP, 

602 

WRITE 

( 

OUP, 

603 


END IF 
END IF 

C 

GO TO 1000 
C 

2000 CONTINUE 


THEN 

) TITLE, NREV, NPAGE 

) 

) 16 , Rl , R2 , GLAT, ELON, HSS , NREC 


C 


IF ( IREC. GT . 3. AND. IPLOT. LE . MAXPLT 

IPLOT * IPLOT +1 
CALL PLOT(X,Y, IREC, RMJD, RNREV) 

END IF 


THEN 


3000 WRITE ( OUP, * ) 
WRITE ( OUP, * ) 
WRITE ( OUP, ★ ) 
WRITE ( OUP, * ) 


NREC, r RECORDS IN FILE. ' 

IPLOT, ' ELEVATION PROFILES PLOTTED, 


CALL UEND 
C 

601 FORMAT ( ' 1 ' , T10 , A20 , 1 OX , ' PASS # 


' , 16 , T90 , 'PAGE ',14,/) 
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602 FORMAT ( ' ' 


JULIAN DAY ' , 2X, ' FRACTION OF DAY',3X, 

ALTIMETER OBS. (M)',5X, ' S/C GEOD LAT 
S/C EAST LONG ',' ELEVATION (M) ',5X, ' RECORD # '/ 

' / 2X, ' ' . 3X, 

nzz z ^ 


' , 5X, ’ — ’ » 

603 rORMAT( ',17 ,T>X, P20 . 18 ~ 3X, F20 . 1775*73 ( 4X, F9 . 6,' 2X) , 5 X,T/ . 3 , 5 x, 110 ) 


STOP 

END 

C 

SUBROUTINE PLOT ( X, Y, IREC, RMJD, RNREV ) 

C 

REAL *4 X(IREC), Y(IREC), RMJD, RNREV 

INTEGER IREC 

CHARACTER* 1 DATE (10) 

C 

CALL USET( 'PERCENTUNITS' ) 

CALL USET ( ' EXTRALARGE ' ) 

CALL UVWPRT (0.0, 99.0, 0.0, 99.0) 

CALL UOUTLN 
C 

C PRINT TITLE, JULIAN DATE AND PASS NUMBER. 

C 

CALL USET ( ' C JUST ' ) 

CALL USET( 'TJUST' ) 

CALL UPSET ( 'PRECISION' ,5.0) 

CALL UPRINT< 50. ,96. , 

6 ' BLACK SEA ALTIMETER ELEVATION PROFILES ' ) 

CALL UPRINT( 50. ,92. , 'MODIFIED JULIAN DATE : S') 

CALL UMOVE(70. ,92. ) 

CALL UPRNT1( RMJD, 'REAL' ) 

CALL UPRINT(50. , 88. , 'PASS : $') 

CALL UMOVE(60. ,88. ) 

CALL UPRNTl (RNREV, 'REAL' ) 

CALL ZTIME( DATE, 8 ) 

CALL FMOVE ( DATE ( 10 ) , 1 , ' $ ' ) 

CALL UPRINT ( 8 8 . , 6 . , DATE ) 

CALL UPRINT ( 88 . , 3 . , ' STX/ZMAYA$ ' ) 

C 

C DRAW AND LABEL AXES (DEFAULT TIC MARKS) 

C 

CALL USET ( 'NOORIGIN' ) 

CALL UVWPRT(5. ,95. , 7. ,90. ) 

CALL USET ( ' DSYMBOL ' ) 

CALL UPSET ( ' SYMBOL ',5.0) 

CALL UPSET (' SZSYMBOL' , 1 . 0) 

CALL USET( 'LARGE' ) 

CALL UPSET ( 'XLABEL' , 'FRACTION OF DAY PAST MIDNIGHTS') 
CALL UPSET ( 'YLABEL' , 'SURFACE ELEVATION IN METERS $ ' ) 
CALL USET( ' XBOTH ' ) 

CALL USET ( ' YBOTH ' ) 

CALL USET ( ' OWNSCALE ' ) 

C CALL UPSET ( 'TICX' ,TCX) 

C CALL UPSET ( ' TICY ' , TCY ) 

XMIN =* X( 1 ) 

XMAX = X ( IREC ) 

YMIN =» Y( 1 ) 

YMAX =* Y < 1 ) 

DO 111 I = 1 , IREC 

YMIN = AMIN 1 ( Y ( I ) , YMIN ) 

YMAX = AMAX1 ( Y(I), YMAX ) 

111 CONTINUE 

SX » 0.05 * ( XMAX - XMIN ) 

SY « 0.10 * ( YMAX - YMIN ) 

XMIN = XMIN - SX 
XMAX = XMAX + SX 
YMIN = YMIN - SY 
YMAX » YMAX + SY 

CALL UWINDO( XMIN, XMAX, YMIN, YMAX) 
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CALL UAXIS ( XMIN, XMAX, YMXN, YMAX) 

C PLOT POLY-LINE 

C 

XPT - FLOAT ( IREC) 

CALL ULINE ( X , Y , XPT ) 

c 

C DISPLAY FOR SCREEN 

C 

c CALL UPAUSE 

C 

c TERMINATED THIS PLOT 

C 

CALL UERASE 
CALL URESET 
C 

RETURN 

END 

//* 

//GO.FT05F001 DO DSN-ZMAYA. BLACK. DATA, DISP-SHR 

//* **** **** .... **.. **** ***. 

//* 

BLACK SEA DATA 

//* **** * * * 

//* 

//GO. FT06F001 DD SYSOUT** 

// EXEC NOTIFYTS 
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noooooononnoooooooonoononoooonnnnnno 


A listing of PROGRAM GRNTRK 


/ /2MAYABAK JOB ( GO 109 , 36 0 , 7 ) , AYAU, TIME- ( 1 , 00 ) , CLASS-O, MSGCLASS-X 
/* JOBPARM LINES-60 
//TEMPLATE EXEC G38PLOT 
//FORT. SYS IN DD * 


FORMAT OF GEOS- 3 


VARIABLE 


TYPE 

14 

12 

12 

14 

14 

14 

R8 

R8 

14 

14 

R4 

14 

14 

14 

R4 

14 

14 

14 

14 

14 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 


/SEAS AT ALTIMETER DATA 
DESCRIPTION 
SATELITE ID 

MEASUREMENT TYPE ( 42* OVER LAND , 43 - OVER WATER ) 
TIME SYSTEM ( NM ) 

STATION NUMBER 
PREPROCESSING INDICATORS 
MODIFIED JULIAN DATE OF OBSERVATION 
FRACTION OF DAY PAST MIDNIGHT (GMT) 

ALTIMETER OBSERVATION 
SATELLITE GEODETIC LATITUDE 
SATELLITE EAST LONGITUDE 
MEASUREMENT STANDARD DEVIATION 
NET INSTRUMENT CORRECTION 
METEOROLOGICAL DATA WORD 
NET MEDIA CORRECTIONS 

GEOID HEIGHT ABOVE REFERENCE ELLIPSOID (METERS) 

NET OCEAN DYNAMIC CORRECTIONS (MM) 

INDICATED SURFACE ELEVATION (MM) 

S/C REVOLUTION NUMBER 

MEAN SEA SURFACE ELEVATION (MARSH/MARTIN '81 (MM)) 
DOD REFERENCE RADIAL ORBIT DIFFERENCE (MM) 

H 1/3 (CM) 

AGC ( DB ) 

WIND SPEED (CM/SEC) 

SURFACE ELEVATION PREPROCESSING WORD 
DRY TROPOSPHERIC CORRECTION (MM) 

FNOC WET TROPOSPHERIC CORRECTION (MM) 

SMMR WET TROPOSPHERIC CORRECTION (MM) 


( METERS ) 

( IE-6 DEGREES) 
( IE-6 DEGREES) 
( METERS ) 

(MM) 

(GEODYN VOL 3) 
(MM) 


IONOSPHERIC CORRECTION 


(MM) 


BAROTROPIC DYNAMIC SEA SURFACE CORRECTION (MM) 


SOLID EARTH TIDE 
SCHWIDERSKI OCEAN TIDE 
PARKE OCEAN TIDE 


(MM) 

(MM) 

(MM) 


C 


C 


C 


c 


PARAMETER ( NMAX * 1000 ) 
PARAMETER ( MAXREC - 10000 ) 


INTEGER* 2 
INTEGER* 4 
INTEGER 
INTEGER 
REAL* 8 
REAL* 4 
REAL *4 
REAL* 4 


12,13,117,118, I 19, 12 0,12 1,122, 123, 124, 125, 126, 127, 120 

II, 14, 15, 16, 17, 18, 19, 110, 111, I 12, I 13, I 14, 115, 116 

IGRID, INP, IREC, ITRK, LPP, NLINE, NPAGE , NREC, OUP 

K75, K78 

R1,R2 

R3,R4 

GLAT, ELON, HSS 

X(NMAX) , Y(NMAX) , XP ( NMAX ) , YP(NMAX) 


INTEGER 
INTEGER 
REAL *4 
REAL *4 
CHARACTER* 1 


NPROJ, LINEH , LINEV 
LIST ( 6 ) 

PENS ( 7) 

PLAT , PLONG , VLAT 1 , VLAT2 , VLONG 1 , VLONG2 
DATE (10) 


DATA PENS / 7*4.0 / 

DATA LIST / 8, 110, 75, 0, 7, 10 / 

DATA NPROJ / 16 / 

DATA PLAT / 90. /, PLONG / 0. / 

DATA VLONG 1 / 26.5 /, VLONG2 / 42.5 / 
DATA VLAT1 / 40. /, VLAT 2 / 48. / 

DATA LINEH / 0 /, LINEV 10 / 

DATA INP / 5 /, OUP / 6 / 
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c 

c. 

c 


DATA 

DATA 


IREC / 0 /, 
IGRID / 1 / 


NREC / 0 /, ITRK / 0 / 


C 

c. 

c 


.PROTECT REGIONAL MAP 

CALL USTART 

CALL UINQES <1.0, SUPRT ) 

2ii s "£; • - r c " 1 1 • 0 ■ LISI ’ 6 • »- io "". ■>»« i 

CALL UVWPRT(0. 0,99. 0,0. 0,99.0) 

CALL UOUTLN 

.PUT ON LABELS 


C 

c. 

c 


CALL USET ( ' C JUST ' ) 

CALL USET( ' T JUST ' ) 

CALL USET( 'YELLOW' ) 

CALL USET ( 'SOFTWARE' ) 

CALL USET ( ' EXTRALARGE ' ) 

CALL UPSET ( 'VERTICAL SIZE ',2.0) 

CALL UPSET ( 'HORIZONT SIZE', 1.6) 

CALL UrONT ( ' TROM ' ) ' 

CALL UPRINT( 50. , 96 . , 

‘iF ( IGRID. G ? E T? /S ™e 5 T GROUND TRACKS OVER BLACK SEA$ '» 
LINEH - IFIX( VLAT2 - VLAT1 ) 
enS INEV - IFIX( VLONG2 - VLONG1 ) 

S CA ^NE^INEV ROJ ' PLAT ' PLONG ' VLAT1 ' ' VLOIWl , VLONG2 , 

•LABEL AXIS AND WRITE LON AND LATS ON PLOT 
IT ( IGRID. GT. 0 ) THEN 


/ LINEV 


420 


C 

C. 

C 


430 


C 

c. 

c 


DLONG * ABS < VLONG2 - VLONG1 ) 

LINES - LINEV + 1 
DIX - 100.0 / LINEV 
YLONGI - -3,5 
DO 420 1*1, LINES , 2 

V » VLONG1 + (DLONG * (I _ i n 
XIaATIT * DIX * (I - 1 ) n 

CALL UMOVE ( XLATIT, YLONGI ) 

CALL UPRNT1(V, 'REAL' ) 

CONTINUE 

.WRITE LABEL 

CALL USET{ 'BJUS ' ) 

CALL UMOVE (50, ,-7.5) 

CALL UPRNTl ( ' EAST LONGITUDE$ ' , ' HORIZ ' ) 

S&VZSffi *?!* VLAT11 ' 

DIY » 100.0 / LINEH 
XLT =-3.5 

DO 430 I ■ 1, LINES, 2 

V = VLAT1 + (DLAT * (I - 1), 

YLO - DIY * (I - 1, ' 

CALL UMOVE ( XLT , YLO ) 

CALL UPRNT1( V, 'REAL' ) 

CONTINUE 

.WRITE LABEL 

CALL USET ( ' MJUS ' ) 

CALL UMOVE (-7.5,50.0) 

CALL UPRNTl ( 'GEODETIC LATITUDES 

END IF 


' VERTI ' ) 
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on o non non non o on o non 


....WRITE DATE AND ID ON PLOT 

CALL ZTIME ( DATE , 8 ) 

CALL FMOVE ( DATE ( 10 ) , 1 , ' $ ' ) 

CALL UPRINT( 79 . , -0 . , DATE) 

CALL UPRINT ( 9 6 . , - 8 . , ' STX/ ZMAYA5 ' ) 

READ ( INP ) 11, 12, 13,14,15, 16, Rl,R2, 17, 18, R3, 

6 19, 110, I11,R4, 112, 113, 114, 115, 116, 117, 118, 119, 

6 120,121,122,123,124,125,126,127,128 

NREV » 114 

IF ( II. EQ. 7502701 ) K75 » K75 + 1 
IF ( II. EQ. 7806401 ) K78 » K78 + 1 
REWIND INP 

DO 111 IJK =* 1, MAX RE C 

READ ( INP, END - 2000 ) I 1 , 12 , 13 , 14 , 15 , 16 , R1 , R2 , 17 , 18 , R3 , 

6 19, I 10, 111, R4, I 12, 113, I 14, 115, I 16, 117, 118, I 19, 

t 12 0,12 1,122, 123, 124, 125, 126, 127, 128 

NREC - NREC + 1 
GLAT - 1.0E-6 * FLOAT (17) 

ELON - 1.0E-6 * FLOAT (18) 

HSS - 1.0E-3 * FLOAT (113) 

IF ( 114. EQ. NREV ) THEN 
IREC - IREC + 1 
X< IREC) » ELON 
.Y ( IREC) » GLAT 
ELSE 

SET UP TO PLOT THE PREVIOUS GROUND TRACK 

IF ( IREC. GT. 1 ) THEN 
ITRK - ITRK + 1 

CALL PLGRTK (X,Y, XP,YP, IREC) 

END IF 
NREV = 114 

IF ( II. EQ. 7502701 ) K75 - K75 + 1 
IF ( II. EQ. 7806401 ) K78 = K78 + 1 
IREC - 1 
X(IREC) = ELON 
Y ( IREC) = GLAT 
END IF 

111 CONTINUE 

000 CONTINUE 

IF ( IREC. GT. 1 ) THEN 

ITRK = ITRK + 1 
CALL PLGRTK ( X , Y , XP , YP , IREC ) 

END IF 

WRITE ( OUP, * ) ' ' 

WRITE ( OUP, * ) ' ', NREC, ' RECORDS IN FILE.’ 

WRITE ( OUP, * ) ' ' 

WRITE ( OUP, * ) ' ITRK, ' GROUND TRACKS PLOTTED, WHERE' 

WRITE ( OUP, * ) ' ' 

WRITE ( OUP, * ) ' ', K75, ' TRACKS FROM GEOS-3 AND' 

WRITE ( OUP, * ) ' ' 

WRITE ( OUP, * ) ' ', K78, ' TRACKS FROM SEASAT.' 

CALL UPAUSE 
CALL UEND 


C 
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SUBROUTINE PLGRTK ( X, Y, XP, YP, IREC ) 

REAL‘4 X(IREC), Y(IREC), XP(IREC), YP(IREC) 
INTEGER IREC 

C 

CALL USET( ' PERCENTUNITS ' ) 

CALL UVWPRT (0.0, 99.0, 0.0, 99.0) 

CALL USET( 'OWNSCALE' ) 

C 

C PLOT GROUND TRACK 

C 

DO 111 I - 1, IREC 

CALL WDPOS ( Y ( I ) , X ( I ) , XP ( I ) , YP ( I ) ) 

111 CONTINUE 
C 

CALL USET( 'LNULL' ) 

XPT » FLOAT (IREC) 

CALL ULINE ( XP , YP , XPT ) 

c 

RETURN 

END 


***★ # i 


//GO.PT05F001 DD DSN-ZMAYA.BLACKA.DATA,DISP-SHR 

//♦ *** + + * + * ttt* #*** **« + # + ** **** +++i 

//* 

BLACK SEA DATA 

//* + *♦* * + + + + * + * ♦ + + * *** + **** 


BLACK SEA DATA 


//GO. FT06F00 1 DD SYSOUT** 

//^XK^oJSfYTS 0 DSN - S *« • W ^ATA2 , DISP-SHR, LABEL- ( , , ,IN) 
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nnnnnnnonoooonnftnnfinonononnnnnnnnonnno 


A listing of PROGRAM XOVERO 


//ZMAYABXO JOB ( GO 109 , 360 , 2 ) , AYAU, TIME- ( 0 , 30 ) , CLASS-O, MSGCLASS-X 
// EXEC TORTVC 
//SYSIN DD * 

FORMAT OF GEOS-3/SEASAT ALTIMETER DATA 


VARIABLE 


# 


TYPE 

DESCRIPTION 




14 

SATELITE ID 




12 

MEASUREMENT TYPE ( 42- OVER LAND 

, 43 ■ 

- OVER 

WATER 

12 

TIME SYSTEM ( NM ) 




14 

STATION NUMBER 




14 

PREPROCESSING INDICATORS 




14 

MODIFIED JULIAN DATE OF OBSERVATION 



R8 

FRACTION OF DAY PAST MIDNIGHT 

(GMT) 



R0 

ALTIMETER OBSERVATION 

( METERS ) 


14 

SATELLITE GEODETIC LATITUDE 

( IE-6 

DEGREES ) 

14 

SATELLITE EAST LONGITUDE 

( IE-6 

DEGREES ) 

R4 

MEASUREMENT STANDARD DEVIATION 

( METERS ) 


14 

NET INSTRUMENT CORRECTION 

(MM) 



14 

METEOROLOGICAL DATA WORD 

(GEODYN VOL 

3) 

14 

NET MEDIA CORRECTIONS 

(MM) 



R4 

GEOID HEIGHT ABOVE REFERENCE ELLIPSOID 

( METERS ) 

14 

NET OCEAN DYNAMIC CORRECTIONS 

(MM) 



14 

INDICATED SURFACE ELEVATION 

(MM) 



14 

S/C REVOLUTION NUMBER 




14 

MEAN SEA SURFACE ELEVATION { MARSH/MARTIN '81 

(MM) ) 

14 

DOD REFERENCE RADIAL ORBIT DIFFERENCE 

(MM) 


12 

H 1/3 

(CM) 



12 

AGC 

(DB) 



12 

WIND SPEED 

(CM/SEC) 


12 

SURFACE ELEVATION PREPROCESSING WORD 



12 

DRY TROPOSPHERIC CORRECTION 

(MM) 



12 

FNOC WET TROPOSPHERIC CORRECTION 

(MM) 



12 

SMMR WET TROPOSPHERIC CORRECTION 

(MM) 



12 

IONOSPHERIC CORRECTION 

(MM) 



12 

BAROTROPIC DYNAMIC SEA SURFACE CORRECTION (MM) 

12 

SOLID EARTH TIDE 

(MM) 



12 

SCHWIDERSKI OCEAN TIDE 

(MM) 



12 

PARKE OCEAN TIDE 

(MM) 




C 


c 


c 


PARAMETER ( MPASS -2 00 ) 
PARAMETER ( NMAX = 200 ) 
PARAMETER ( SIGE * 2.0) 


INTEGER* 2 
INTEGER* 4 
REAL* 8 
REAL* 4 
INTEGER* 4 
REAL *4 
REAL *4 
REAL* 4 
REAL* 8 
REAL* 8 
REAL* 8 
REAL *4 
INTEGER* 4 
INTEGER* 4 
INTEGER 


12,13,117,118,119,120,121,122,123,124,125,126,127,128 
II, 14, 15, 16, 17, 18, 19, 110, 111, 112, 113, 114, 115, 116 


R1,R2 

R3,R4 

IDSAT 

GLAT , ELON 

X(NMAX) , Y(NMAX) , Z(NMAX), YR(NMAX) 

GX ( NMAX ) , GY ( NMAX ) , DDZ ( NMAX ) , D(NMAX) 

F (NMAX, 3), FT (3, NMAX) 

G < NMAX ) 

A( MPASS ) , B ( MPASS ) , C(MPASS) 

SLON 1( MPASS ) , SLAT 1( MPASS ) , SLON2 ( MPASS ) , SLAT2 ( MPASS ) 
IDS (MPASS), NREV( MPASS), NBAD( MPASS) 

ICX ( MPASS, MP ASS ) , ISUM(MPASS), LN ( MPASS ) , LY ( MPASS ) 

IDP , INP, IREC, NPASS , NREC , OUP 


DATA INP / 3 /, OUP / 6 / 

DATA IREC / 0 /, NREC / 0 /, NPASS / 0 / 


READ ( INP 

& 

& 


II, 12, 13, 14, 15, I6,R1,R2, 17, I8,R3, 

19, 110, 111, R4, I 12, 113, I 14, 115, 116, I 17, 118, I 19, 
120, 121, 122, 123, 124, 125, 126, 127, 128 
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non 


IDP - 114 
REWIND INP 
C 

WRITE ( OUP, CURVE FITTING FOR EACH SATELLITE PASS ' 

WRITE ( OUP, * ) ' ' 

C 

1000 READ ( INP, END - 2000 ) 11,12,13,14, 15, 16, R1,R2, 17,18, R3, 

* 19, 110, II 1,R4, 112, 113, I 14, 115 ,116, 117, 118, 119, 

* 120,121,122,123,124, 125,126 , 127 , 128 

IT ( 114. EQ. 12793 ) GO TO 999 


NREC - NREC + 1 
GLAT » 1.0E-6 * FLOAT (17) 
ELON ■ 1.0E-6 * FLOAT (18) 
HSS » 1.0E-3 * FLOAT (113) 
IDSAT - II 

IF ( 114. EQ. IDP ) THEN 
IREC - IREC + 1 
IDP - 114 
X(IREC) - ELON 
Y(IREC) - GLAT 
2 (IREC) - HSS 
ELSE 


.DETERMINE THE SECOND-ORDER POLYNOMIAL FOR A PASS 


IF ( IREC. GE. 5 ) 


THEN 


NPASS - NPASS + 1 
NREV(NPASS) - IDP 
IDS (NPASS) » IDSAT 
SLONl( NPASS) 

SLAT1( NPASS) 

SLON2 ( NPASS ) 

SLAT 2 ( NPASS ) 


CALL CURVE ( 


C 

C 


X( 1) 

Y(l) 

* X(IREC) 

- Y(IREC) 

NPASS, X, Y, Z, YR, IREC, NREV(NPASS), 

IDS ( NPASS ) , A ( NPASS ) , B(NPASS), C(NPASS), 

NBAD ( NPASS ) , F, FT, G, GX, GY, DDZ , D, SIGE ) 


999 


END IF 

IDP - 114 
IREC - 1 
X(IREC) - ELON 
Y ( IREC ) ■ GLAT 
Z(IREC) ■ HSS 
END IF 

CONTINUE 
GO TO 1000 


2000 CONTINUE 

IF ( IREC. 


GE. 5 ) 


THEN 


NPASS - NPASS + 1 
NREV( NPASS) * IDP 


IDS (NPASS) 
SLONl( NPASS) 
SLAT1( NPASS) 
SLON2 ( NPASS ) 
SLAT2 ( NPASS ) 
CALL CURVE ( 


C 

c 


IDSAT 
- X(l) 

» 7(1) 

= X ( IREC ) 

= Y(IREC) 
NPASS, X, Y, 
IDS (NPASS) , 
NBAD (NPASS) , 


Z, YR, IREC, NREV( NPASS), 

A( NPASS ) , B ( NPASS ) , C(NPASS), 

F, FT, G, GX, GY, DDZ, D, SIGE) 


END IF 
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n non 


3000 WRITE< OUP 
WRITE! OUP 
WRITE! OUP 
WRITE ( OUP 
C 

CALL XOVER 

6 

C 

STOP 

END 

C 

SUBROUTINE 

& 

C 

REALM 
REAL* 4 
REAL* 8 
REAL *8 
REAL* 8 
REAL* 8 
INTEGER* 4 
INTEGER* 4 
INTEGER* 4 


/ * ) ' ' 

, * ) ' NREC, ' RECORDS IN TILE.' 

r * )' ' 

, * ) ' NPASS, ' PASSES REPRESENTED BY AX**2+BX+C 

( NPASS, NREV, IDS, A, B, C, SLON1, SLATl, 

SLON2, SLAT 2 , ICX, ISUM, LN, LY, NBAD ) 


CURVE ( NPASS, X, Y, Z, YR, IREC, IDP, IDSAT, A, B, 
NB, T, FT, G, GX, GY, DDZ, D, SIGE ) 

X(IREC), Y(IREC), Z(IREC), YR(IREC) 

GX(IREC), GY (IREC), DDZ (IREC), D(IREC) 

F ( IREC , 3 ) , FT( 3 , IREC ) 

E( 3 , 3 ) , El ( 3 , 3 ) , AE ( 6 ) 

G(IREC), H( 3 ) 

A, B, C, DET, YD 

IDP, IDSAT, IREC, NB, NPASS 

IDIAG<3) / 1, 3, 6 / 

OUP / 6 / 


C, 


EVICT THE BAD DATA POINTS BY COMPARING THE ABSOLUTE DIFFERENCE 
BETWEEN ADJACENT DATA POINTS 


NB * 0 

IREC 1 » IREC - 1 
C 

DO 111 I - 1, IREC 
GX ( I ) » -1.0 

111 CONTINUE 
C 

DO 112 I - 1, IREC 1 

DDZ ( I ) = ABS ( Z(I+1) - Z(I) ) 

112 CONTINUE 
C 

DO 113 J ■ 1, IREC 1 

IF(DDZ( J) . LT.SIGE) THEN 
K = J 
GO TO 200 
END IF 

113 CONTINUE 
C 

200 CONTINUE 

DO 222 J = 1, IREC 

IF(DDZ( IREC -J) .LT.SIGE) THEN 
IRECG * IREC-J 
GO TO 300 
END IF 
222 CONTINUE 
C 

300 KB = 0 

DO 333 1 = 1, IRECG- 1 

IF( DDZ ( I ) .GT. SIGE) KB = KB + 1 
333 CONTINUE 
C 

LG = 1 

GX ( LG ) = X ( K ) 

GY ( LG ) = Y ( K ) 

C 

400 LG = LG + 1 
K = K + 1 

IF ( DDZ ( K-l ) .GT. SIGE) THEN 
K = K + 1 

DO 444 1=1, KB- 1 

IF ( DDZ ( K- 1 ) . GT . SIGE ) K = K + 1 
444 CONTINUE 
END IF 
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GX(LG) » X< K ) 

GY(LG) - Y(K) 

IT ( K . LE . IRECG ) GO TO 400 

DO 555 I - 1, LG 

F(I, 1) « GX(I)*GX(I) 
F(I,2) - GX( I ) 

F( I, 3) - 1.0D0 
G(I) - GY ( I ) 

555 CONTINUE 


C 


c 


DO 666 I - 
DO 667 J 
FT(J,I) 

667 CONTINUE 
666 CONTINUE 


1, LG 

- 1, 3 

- F<I,J> 


DO 777 I - 
DO 778 J 
E(I,J) - 
DO 779 

779 CONTINUE 

778 CONTINUE 
777 CONTINUE 


1, 3 
- 1, 3 
' 0 . DO 
K - 1, LG 
- E(I,J) 


+ 


FT(I,K)*F(K, J) 


C 


c 


DO 888 I ■ 1, 3 

H( I ) - 0 . DO 
DO 889 K » 1, LG 
„„„ H(I) - H( I ) + FT (IjK) *G( K 

889 CONTINUE 
888 CONTINUE 


AE( 1) * E(l, 1) 
AE(2) - E( 1,2) 
AE(3) - E( 2 , 2 ) 
AE(4) ■ . E( 1 , 3 ) 
AE(5) - E( 2 , 3 ) 
AE(6) - E( 3, 3 ) 


C 


CALL DSOLVE ( AE, H, IDIAG, 3, 


. TRUE . , 


. TRUE . ) 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


A - H( 1 ) 
B - H( 2 ) 
C - H( 3) 


DET ■ E ( 1 , 1 ) * ( 
4 E( 1 , 2 ) * ( 
4 E( 1,3) * ( 


E( 2 , 2 ) +E( 3, 3 ) - 
E ( 1,3) *E (2,3) - 
E< 1, 2 ) *E( 2 , 3 ) - 


E( 2, 3 ) *E( 2,3) ) + 
E ( 1 f 2 ) *E( 3,3) ) + 
E( 1,3)*E<2,2) ) 


El (1,1) 
El ( 2 , 2 ) 
£1(3,3) 
EI( 2 , 3 ) 
BI(1. 3) 
EI< 1,2) 


( E(2,2)*E(3,3) 
( E(l, 1)*E(3,3) 
( E( 1, 1)*E<2,2) 
( E( 1 , 2 ) *E( 1 , 3 ) 
( E( 1,2) *E( 2,3) 
( E( 1 , 3 ) *E ( 2 , 3 ) 


E(2,3)*E(2,3) ) / DET 
E( 1 , 3 ) *E( 1 , 3 ) ) / DET 
E(1,2)*E(1,2) ) / DET 
E( 1# 1 ) *E( 2, 3 ) ) / DET 
E( 1 , 3 ) *E( 2 , 2 ) ) / DET 
E( 1,2)*E(3,3) ) / DET 


A 

B 

C 


( El ( 1 , 1 ) *H( 1 ) 
( El ( 1 , 2 ) *H( 1 ) 
( El ( 1, 3 ) *H( 1 ) 


♦ EI< 1 , 2 ) *H ( 2 ) 
+ El ( 2 , 2 ) *H( 2 ) 
+ El ( 2 , 3 ) *H( 2 ) 


+ EI( 1 , 3 ) *H( 3 ) ) 
+ EI( 2 , 3 ) *H( 3 ) ) 
+ El ( 3 , 3 ) *H ( 3 ) ) 


RMS =■ 0. 

DO 999 I = 1, LG 

YD = A*GX(I)*GX(I) + B*GX ( I ) + C 
YR< I ) = SNGL( YD) 

D( I ) “ YR( I ) - GY ( I ) 

RMS * RMS + D(I)*D(I) 

999 CONTINUE 

RMS = SQRT ( RMS /FLOAT (LG) ) 
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n n 


IF ( RMS. GT. 1.00 ) NB - 1 
C 

WRITE ( OUP, 601 ) NPASS, IDP, IDSAT, A, B, C, 

* GX ( 1 ) , GY(1), GX(LG), GY(LG), RMS, LG, NB 

IF ( RMS. GT. 0.005 ) THEN 
WRITE ( 6, * ) 

WRITE ( 6, 602 ) (M,GX(M) ,GY(M) , YR(M) ,D(M) ,M-1,LG) 

WRITE ( 6, * ) ' ' 

END IF 
C 

601 FORMAT ( 1H0, //, IX, 13 , 17 , 19 , 3E14 . 5 , 5F10 . 4 , 16 , 13 ) 

602 rORMAT((lH , T5 , 14 , 3X, 3F10 . 4 , 2X, F10 . 4 ) ) 

C 

RETURN 

END 

C 

SUBROUTINE XOVER ( NPASS, NREV, IDSAT, A, B, C, SLONl, SLAT1, 

6 SLON2, SLAT2, ICX, ISUM, LN, LY, NBAD ) 

C 

INTEGER NREV( NPASS), IDSAT{ NPASS ) , NBAD( NPASS) 

INTEGER ICX (NPASS, NPASS ) , ISUM (NPASS), LN( NPASS), LY( NPASS) 
INTEGER INDEX ( 2 ) 

REALM SLAT1 ( NPASS ) , SLAT2 ( NPASS ) , SLONl ( NPASS ) , SLON2 ( NPASS) 

REALM A( NPASS ) , B(NPASS), C(NPASS) 

REALM T( 2 ) 

REAL *8 DA, DB, DC, RT, SRT, P, Q 

INTEGER ITOTAL, IXM, MAXIX, MX, MXS, MX75, MX78, NPASS, OUP 

DATA INP / 5 /, OUP / 6 / 

DATA MAXIX / 0 / , MX / 0 / , MX75 / 0 /, MX78 / 0 / 

DATA ITOTAL / 0 / 

C 

DO 123 I ■ 1, NPASS 
ISUM(I) - 0 
DO 234 J - 1, NPASS 
ICX(I,J) * 0 
234 CONTINUE 
123 CONTINUE 
C 

WRITE ( OUP, 600 ) 

C 

DO 222 1*1, NPASS- 1 

C 

ir ( NBAD(I). EQ. 0 ) THEN 
C 

WRITE ( OUP, 601 ) NREV ( I ) , IDSAT(I), SLONl(I), SLATl(I), 

A SLON2 ( I ) , SLAT2 ( I ) 

C 

DO 333 J * 1+1, NPASS 
C 

IF ( NBAD(J). EQ. 0 ) THEN 
C 

INDEX ( 1 ) * 0 
INDEX ( 2 ) =* 0 
DB - B( J) - B( I ) 

DA = A( I ) - A( J) 

DC * C ( I ) - C ( J ) 

RT = DB*DB - 4 . 0 * DA* DC 
C 

IF ( RT. GE. 0.0 ) THEN 
IF ( RT. EQ. 0.0 ) THEN 
SRT = 0.0 
ELSE 

SRT = DSQRT(RT) 

END IF 

T ( 1 ) = (DB+SRT)/(2.0*DA) 

T( 2 ) = ( DB-SRT) / ( 2 . 0*DA) 

WRITE ( 6 , * ) T, RT, SRT 

DO 135 L = 1, 2 

IF ( SLON2 ( I ) . GT. SLON 1(1) ) THEN 
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4 


4 


& 


« 


135 

C 


4 

4 


246 

C 


IF( T(L). GE. SLON 1(1). AND. T(L). 
INDEX(L) - INDEX ( L ) + 1 

ELSE 

IF( T(L). LE. SLON 1(1). AND. T(L) . 
INDEX(L) - INDEX ( L) + 1 
END IF 

ir ( SLON2 ( J) . GT. SLON 1 ( J ) ) THEN 
IF ( T(L). GE. SLON 1 ( J ) . AND. T(L) . 
INDEX ( L ) * INDEX ( L ) + 1 

ELSE 

IF ( T(L) . LE. SLONl(J). AND. T(L). 
INDEX(L) - INDEX(L) + 1 
END IF 
CONTINUE 


LE . SLON 2 ( I ) ) 

GE. SLON2 ( I ) ) 


LE. SLON2 ( J) ) 
GE . SLON 2 ( J ) ) 


DO 246 L = 1, 2 
IF ( INDEX(L). 
MX ■ MX + 1 
P * T( L ) 

Q * A(I)*P*P 
WRITE ( OUP, 


EQ. 2 ) THEN 


+ B ( I ) *P + C(I) 

602) NREV(J), IDSAT(J) 
SLON1 ( J) , SLAT1 ( J) , 
SLON2 ( J ) , SLAT2 ( J ) 
ICX(I,J) - ICX( I, J) +1 
ICX(J,I) - ICX(I,J) 

IF ( IDSAT(I). EQ. IDSAT(J) ) 


IF ( 
IF ( 
END IF 
END IF 
CONTINUE 


IDSAT( I) 
IDSAT ( I ) 


EQ. 

EQ 


7502701 

7006401 


THEN 
MX75 
MX7 8 


# P # Q, 


- MX75 + 1 
= MX 7 8 + 1 


C 


END IF 


END IF 
C 

333 CONTINUE 

C 


END IF 

C 

222 CONTINUE 
C 

WRITE ( OUP, * ) 
WRITE ( OUP, * ) 
WRITE ( OUP, * ) 
WRITE ( OUP, * ) 
WRITE ( OUP, * ) 
WRITE ( OUP, * ) 
WRITE ( OUP, * ) 
MXS - MX - MX75 - 
WRITE ( OUP, * ) 
WRITE ( OUP, * ) 


' TOTAL # OF POSSIBLE CROSSOVERS * ' , MX 

TOTAL 4 OF GE03-GE03 CROSSOVERS = ', MX75 

' TOTAL # OF SEASAT-SEASAT CROSSOVERS = ’, MX78 
MX7 8 

0 9 

TOTAL # OF GEOS3-SEASAT CROSSOVERS =* ' , MXS 


456 

345 


C 


DO 345 1 = 1, NPASS 

DO 456 J = 1 , NPASS 

ISUM(I) = ISUM< I ) + ICX ( X , J ) 
CONTINUE 

ITOTAL = ITOTAL + ISUM(I) 
CONTINUE 

ITOTAL = ITOTAL/2 


DO 567 1=1, NPASS 

MAX IX = MAXO ( MAX IX , 
567 CONTINUE 
C 


ISUM(I) ) 


DO 678 1=1, NPASS 

IF ( ISUM(I). EQ. MAXIX 
678 CONTINUE 
C 


IXM = I 


WRITE ( OUP, 603 ) 
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nnaon 


c 


c 


WRITE ( OUP, 604 ) (I, NREV) I), I - 1, NPASS ) 

WRITE ( OUP, 605 ) NPASS 

WRITE ( OUP, 606 ) ( NREV)I),<ICX(I,J),J*1,94),ISUM(I), 

* 1*1, NPASS ) 

WRITE ( OUP, 607 ) NREV(IXM), HAXIX, ITOTAL 

MM - 0 
NN ■ 0 

DO 789 I - 1, NPASS 

ir ( ICX( IXM, I ) . EQ. 1 ) THEN 
MM - MM + 1 
LY(MM) » I 
ELSE 

NN = NN + 1 
LN ( NN ) - I 
END IF 
789 CONTINUE 


WRITE ( OUP, 608 ) 
WRITE ( OUP, 609 ) 


NREV(IXM), MM, NN, MM+NN 
( NREV(LY) I) ) , 1*1, MM ) 


DO 987 J » 1, NN 

WRITE ( OUP, 610 ) NREV ( LN ( J ) ) 

DO 876 I - 1, NPASS 

IF ( ICX( I , LN( J) ) . EQ. 1 ) THEN 
WRITE ( OUP, 611 ) NREV ( I ) 

GO TO 888 
END IF 

876 CONTINUE 

888 CONTINUE 
987 CONTINUE 


WRITE ( 6, * ) ' ' 

WRITE ( 6, 612 ) 

DO 999 1*1, NPASS 

DO 998 J * 1, NPASS 
IF ( ICX(I,J). GT. 

998 CONTINUE 

999 CONTINUE 


1 ) WRITE ( 6, * ) NREV ( I ) , NREV ( J ) 


501 

600 

601 

602 

603 

604 

605 

606 

607 

608 

609 

610 
611 
612 


FORMAT (I5,I10,2F8.3,4F10.6) 

FORMAT) 1H1,T20, 'BLACK SEA CROSSOVER LISTING:',//) 

FORMAT) 1H0,//,T5, 'PASS * ',I6,5X, 'SAT ID * ' ,110, 5X, 4F10.6//) 
FORMAT) 1H , T10 , ' XPASS * ',I6,5X,'SAT ID * ' , 110, 5X, 2F10. 3, 

<. 5X, 4F10 . 6 ) 

FORMAT) 1H1 , T20, ' ROW-COLUMN ELEMENTS OF PASS CORRELATION MATRIX; 

* //) 

FORMAT) ( 1H ,10(3X, 14,16) ) ) 

FORMAT) 1H1,T20, 'PARTIAL (94 COLUMNS) PASS CORRELATION MATRIX ' 
i 'FOR ',15,' PASSES:',//) 

FORMAT ( ( 1H ,10X,I8,5X,94I1,5X,I3) ) 

FORMAT) 1H0, //,T10, 'PASS # ',16,' HAS THE MAXIUM OF ',15, 
i ' CROSSOVERS. ', //,T10, 'THE TOTAL # OF CROSSOVERS IS:’, I* 

FORMAT) 1H1,T5, 'REFERENCE PASS = ' , 16 , T1 10 , 316 , // , 

*■ T10, 'THE DIRECT CROSS-PASSES ARE:' //) 

FORMAT) ( 1H0, 2016 ) ) 

FORMAT) 1H0, //,T5, 'PASS: ',16,' IS INDIRECTLY RELATED TO THE') 

FORMAT ( 1H ,T10, 'REFERENCE PASS VIA PASS ',16) 

FORMAT) 1H1 , T10 ,' POSSIBLE COMBINATIONS OF PASSES THAT YIELD ' 

4 'BAD CROSSOVERS:’,//) 


RETURN 

END 


SUBROUTINE DSOLVE ( A, B , IDIAG, NEQ, FACT, BACK) 


Compute the U*-»T * D * U factorization of the symmetric matrix 
stored m A, if FACT = TRUE; and solve A * X = B if BACK = TRUE. 
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uuuuuuuuuuuuuuu uuu 


A 


B 


ID I AG 
NEO 

FACT 

BACK 


Contains the compacted-col limn form of the upper trianaular 

Right-hand-aid* vector. After backsubetitution, it 
contains the solution. 


Addresses of the diagonal terms in A. 
Number of equations 


If FACT ■ TRUE, then factor A; otherwise do 
If BACK * TRUE , reduce B and backsubstitute: 
do not solve the equations. 


not factor A. 
otherwise 


IMPLICIT REAL *8 

LOGICAL 

DIMENSION 


<A-H,0-Z) 

FACT, BACK 

A(1),B(1), IDIAG ( 1) 


Factor A, reduce B 
JR » 0 . 

DO 400 J » 1, NEQ 
JD - IDIAG(J) 

JH * JD - JR 
IS « J - JH + 2 
C 

IF (JH .LT. 2) GOTO 390 
C 

IF (FACT) THEN 
C 

IF (JH .GT. 2) THEN 

C.. Reduce^coluran J rows is to J-l: do not divide by row diagonal 

ID - IDIAG ( IS - 1) 


100 

C 

c. . 
c 


DO 100 I * IS, J-l 
IR - ID 
ID « IDIAG(I) 

IH - MIN ( ID-IR-1, I-IS+1) 

K F - ( K H + G r 0) A<K) = A<K> ~ *OT(A<K-IH),A<ID-IH>,IH) 

CONTINUE 

ENDIF 

Divide by row diagonal, and reduce diagonal term in column J 


200 

C 


IR * JR + 1 
K * J - JD 

DO 200 I = IR, JD- 1 
ID * IDIAG( K+I ) 

IF ( A< ID) . EQ. 0 . 0 ) GOTO 200 
D - -A( I) 

A( I ) - A( I ) /A( ID) 

A( JD) - A( JD) + D*A( I) 
CONTINUE 


ENDIF 

C 

C.. Reduce RHS 

c IF < BACK > B(J) = B(J) - DOT(A(JR+l),B(IS-l),JH-l) 

390 JR = JD 
400 CONTINUE 
C 
C 


IF (.NOT. BACK) RETURN 
C 
C 

C Divide by diagonal pivots 
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c 

DO 700 I - 1 , NCQ 
ID * IDIAG(I) 

IF <A< ID) .NE.0.0) B( I) 
700 CONTINUE 
C 

C Backsubstitute 

c 

J - NEQ 

JD - IDIAG(J) 

C 

801 D - -B{ J ) 

J - J - 1 

IF (J.LE.O) RETURN 
C 

JR * IDIAG(J) 

IF (JD-JR.GT.l) THEN 
IS * J - JD + JR + 2 
K * JR - IS + 1 
DO 810 I * IS, J 
810 B ( I ) - B { I ) + A( I+K ) *D 
ENDIF 
C 

JD - JR 
GOTO 801 
C 


B ( I ) /A( ID ) 


END 

C 

FUNCTION DOT ( A, B , N ) 

C 

C Compute the dot product of the two N-vectors A and B. 

C 

INTEGER N, I 

REAL* 8 DOT, A(1),B<1) 

C 

DOT*0 . 0 
C 

DO 100 1*1, N 

100 DOT=DOT + A(I)*B<I) 

C 

RETURN 

END 

//* 

// EXEC LINKGOV, REGION*3000K 
//♦SYSLIB DD DSN-SYS2 . IMSLD , DISP*SHR 

//* DD DSN*SYS2 . IMSLS, DISP*SHR 

//* 

//GO.FT03F001 DD DSN=Z MAYA . BLACK . DATA, D I SP*S HR 

//* **** **** **** ***♦ ***♦ ♦*** **** 

//* 

//* BLACK SEA DATA 

//* 

//* * + + + + + ** + ** + ** + * * + + * + + + * + * + * 

//* 

//GO. FT06F00 1 DD SYSOUT=* 

// EXEC NOTIFYTS 
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A listing of PROGRAM XOVER 


c 

c 

c 

c 

c 

c 

c 

C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
. c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


//ZMAYABXW JOB ( GO 109 , 360, 2 ) , AYAU,TIME*( 1,30), CLASS-O, MSGCLASS-X 
/♦JOBPARM LINES-60 
// EXEC FORTVC 
//SYSIN DD * 

C 

FORMAT OF GEOS- 3 /SEAS AT ALTIMETER DATA 


VARIABLE 


TYPE 

14 

12 

12 

14 

14 

14 

R8 

R8 

14 

14 

R4 

14 

14 

14 

R4 

14 

14 

14 

14 

14 

12 

12 

12 

12 

12 

12. 

12 

12 

12 

12 

12 

12 


DESCRIPTION 

SATELITE ID 
MEASUREMENT TYPE ( 

TIME SYSTEM ( NM ) 

STATION NUMBER 
PREPROCESSING INDICATORS 
MODI TIED JULIAN DATE OF OBSERVATION 


42* OVER LAND , 43 * OVER WATER 


(GMT) 

( METERS ) 

( IE-6 DEGREES) 
(IE-6 DEGREES) 
( METERS ) 

(MM) 

(GEODYN VOL 3) 
(MM) 

( METERS ) 


FRACTION OF DAY PAST MIDNIGHT 
ALTIMETER OBSERVATION 
SATELLITE GEODETIC LATITUDE 
SATELLITE EAST LONGITUDE 
MEASUREMENT STANDARD DEVIATION 
NET INSTRUMENT CORRECTION 
METEOROLOGICAL DATA WORD 
NET MEDIA CORRECTIONS , 

GEOID HEIGHT ABOVE REFERENCE ELLIPSOID 
NET OCEAN DYNAMIC CORRECTIONS (MM) 

INDICATED SURFACE ELEVATION (MM) 

S/C REVOLUTION NUMBER 

MEAN SEA SURFACE ELEVATION (MARSH/MARTIN '81 (MM)) 
DOD REFERENCE RADIAL ORBIT DIFFERENCE (MM) 

H 1/3 (CM) 

AGC ( DB ) 

WIND SPEED (CM/SEC) 

SURFACE ELEVATION PREPROCESSING WORD 
DRY TROPOSPHERIC CORRECTION (MM) 

FNOC WET TROPOSPHERIC CORRECTION (MM) 

SMMR WET TROPOSPHERIC CORRECTION (MM) 

IONOSPHERIC CORRECTION (MM) 

BAROTROPIC DYNAMIC SEA SURFACE CORRECTION (MM) 
SOLID EARTH TIDE (MM) 

SCHWIDERSKI OCEAN TIDE (MM) 

PARKE OCEAN TIDE MM) 


PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

INTEGER* 2 
INTEGER* 4 
REAL *8 
REAL *4 
INTEGER* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 


SIGMAG 
SIGMAS 
XDMAX 
ETOP * 
EBOT * 
SIGE = 
MAXOVR 
MPASS 
NMAX 
MS YM 


0.25 
0.10 
100 . 

100 . 

0. 

2 . 

2100 
142 
140 

{ MPASS * ( MPASS + 1 ) 


/2 ) 


12,13,117,118,119,120,121,122,123,124,125,126,127,128 

I 1, I 4, I 5, I 6, I7 ,I8,I9, I1 °,IH,I12,I13,I14,I15,I16 

R1,R2 

R3,R4 

IDSAT 

GLAT, ELON, HSS 

XI ( NMAX) , YI ( NMAX ) , 21 ( NMAX ) , 2G( NMAX) , YR( NMAX) 

SX(NMAX) , SY ( NMAX) , SZ ( NMAX ) , SZG( NMAX) , DDZ ( NMAX) , D( NMAX) 
AX (MAXOVR, MPASS ) , BX(MAXOVR), XVR ( MPASS , MPASS ) 

W( MAXOVR), AXT( MPASS, MAXOVR) , SCR ( MPASS , MPASS ) 

AXTBX ( MPASS ) , AXTAX(MSYM) 

X (MPASS, NMAX ) , H ( MPASS , NMAX) , HG( MPASS , NMAX ) 

SLON 1 ( MPASS ) , SLAT 1 ( MPASS ) , SLON2 ( MPASS ) , SLAT 2 ( MPASS ) 
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non 


c 


c 


c 


REAL* 4 
REAL* 4 
REAL* 8 
REAL* 8 
REAL* 8 
REAL* 8 
INTEGER* 4 
INTEGER* 4 
INTEGER* 4 
INTEGER 
INTEGER 


U(MPA SS), S ( HP ASS ) , CIMSL< NMAX-1 , 3 ) 

EX ( HP ASS , MPASS ) , GY ( MPASS , MPASS ) 

DAXTBX ( MPASS ) , DAXTAX ( MSYM ) , ERCOVM ( MS YM ) 

F ( NMAX , 3 ) , FT ( 3 , NMAX ) 

G(NMAX) 

A( MPASS ) , B ( MPASS ) , C(MPASS) 

IDS (MPASS), IDIAG( MPASS ) , NREV( MPASS), NIREC( MPASS) 
ICB( MPASS), ISUM( MPASS), LN( MPASS), LY( MPASS) 

ICX( MPASS, MPASS) 

IDP, IREC, IRECA, IXM, NPASS, NREC 
ING, INP, OUP, OUG, OUC 


DATA INP / 3 /, ING / 2 /, OUG / 8 /, OUP / 6 /, OUC / 9 / 

DATA IREC / 0 /, NREC / 0 /, NPASS / 1 / 


G2 » SIGMAG*SIGMAG 
S2 « S IGMAS * SIGMAS 
SMGG - 1 . 0/ ( 2 . *G2 ) 

SMSS - 1.0/(2.*S2) 

SMSG - 1.0/(S2+G2) 

READ ( INP ) 11,12,13,14,15, I6,R1,R2, 17, 18, R3, 

6 19,110, 111, R4, 112, 113, 114, 115, 116, 117, 118, 119, 

& 120,121,122,123, 124,125, 126,127, 128 

IDP » 114 
REWIND INP 


C 

1000 READ 
& 

& 

C 

READ 

& 

C***** 

IF ( 

C***** 


( INP, END * 2000 ) II, 12, 13, 14, 15, 16, R1,R2, 17, 18, R3, 

19,110, 111, R4, 112, 113,114, 115, 116,117, 118, 119, 
120,121,122,123,124,125,126,127,128 

( ING, 501, END - 2000 ) IDUMl, GL, EL, DUM2 , DUM3 , DUM4 , 

DUM5 , DUM6 , GUND 

114. EQ. 12793 ) GO TO 999 


NREC » NREC + 1 
GLAT * 1.0E-6 * FLOAT ( 17 ) 
ELON - 1.0E-6 ♦ FLOAT (18) 
HSS - 1.0E-3 * FLOAT (113) 
IDSAT -II 


C 


IF ( 114. EQ. IDP ) THEN 

IF ( HSS. GT. EBOT. AND. HSS. LT. ETOP ) THEN 

IREC - IREC + 1 
IDP - 114 

X( NPASS, IREC) =» ELON 
H( NPASS, IREC) - HSS 
XI (IREC) - ELON 
YI(IREC) =* GLAT 
ZI(IREC) = HSS 
ZG(IREC) = GUND 
END IF 


C 


ELSE 


DETERMINE THE SECOND-ORDER POLYNOMIAL FOR A PASS 


C 


IF ( IREC. GE. 5 ) THEN 


& 

& 


NREV( NPASS) 
IDS (NPASS) = 
CALL CURVE ( 


NIREC( NPASS) 
SLONl( NPASS) 
SLAT1( NPASS) 
SLON2 ( NPASS) 


= IDP 
IDSAT 

NPASS, XI, YI, ZI 
IDSAT, A( NPASS), 
F, FT, G, SX, SY, 
= IREC 
= XI(1) 

* YI(1) 

= XI (IREC) 


, ZG, YR, 
B( NPASS) , 
SZ, SZG, 


IREC, IDP, 

C( NPASS) , 
DDZ, D, SIGE 


) 
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100 

c 


SLAT2 ( NPASS ) 
DO 100 M - 1 
X { NPASS , M ) 

H ( NPASS , M ) 
HG( NPASS, M) 
CONTINUE 


YI ( IREC ) 
IREC 
XI(M) 
ZI(M) 
ZG(M) 


C 


ELSE 


NPASS 


NPASS - 1 


END IF 

C 

IDP - 114 
IREC - 1 

NPASS - NPASS + 1 
IF ( HSS. GT. EBOT. 

X( NPASS, IREC) » ELON 
H( NPASS, IREC) - HSS 
XI (IREC) - ELON 
YI(IREC) - GLAT 
ZI(IREC) m HSS 
ZG(IREC) ■ GUND 
ELSE 

IREC - IREC - 1 
END IF 
END IF 

C 

999 CONTINUE 
GO TO 1000 

C 

^2000 CONTINUE 

^ IF ( IREC, GE , 5 ) THEN 


AND . HSS . LT , ETOP ) THEN 


& 

6 


NREV( NPASS) 
IDS (NPASS) ~ 
CALL CURVE ( 


101 


C 

C 

C 


NIREC( NPASS) 
SLONl( NPASS) 
SLAT1( NPASS) 
SLON2 ( NPASS ) 
SLAT2( NPASS) 
DO 101 M - 1 
X ( NPASS , M) 
H( NPASS, M) ■ 
HG( NPASS, M) 
CONTINUE 

ELSE 

NPASS = NPASS - 
END IF 


> IDP 
IDS AT 
NPASS , 
IDSAT, 

F, FT, 

■ IREC 

- XI(1) 

- YI ( 1 ) 

■ XI (IREC) 

- YI(IREC) 
, IREC 

■ XI ( M ) 

■ ZI ( M) 

* ZG(M) 


XI, YI, SI, ZG, YR, 
A( NPASS), B( NPASS), 
SY, SZ, SZG, 


G, SX, 


IREC, IDP, 

C ( NPASS ) , 

DDZ , D, SIGE 


3000 WRITE( OUP, * ) 
WRITE ( OUP, * ) 
WRITE ( OUP, * ) 
WRITE ( OUP, * ) 


NREC, ' RECORDS IN FILE.’ 

NPASS, ' PASSES REPRESENTED BY AX**2+BX+C' 
CALL XOVER , »REV IDS A B, C, SLON1 . SLAT1 , 


CT ° r SbUNI, SLAT 1 , 

SLON2 , SLAT2 , ICX, EX, GY, ISUM, LN , LY, IXM ) 


CALL CLR ( XVR, MPASS*MPASS) 
CALL CLR ( AX, MAXOVR*MPASS ) 
CALL CLR ( BX, MAXOVR) 

CALL CLR ( AXTAX , MSYM ) 

CALL CLR ( AXTBX, MPASS ) 
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c 

c 

201 

202 


203 

& 

6 


C 

c & 

204 


& 

205 


CALL DCLR { DAXTAX, MSYM) 

CALL DCLR ( DAXTBX r MPASS ) 

WRITE { OUP, 600 ) 

WRITE ( OUP, START INTERPOLATION FOR EACH PASS' 

DO 200 I - 1, NPASS 
IREC - NIREC(I) 

IF ( X(I,1). LT . X<I,IREC) ) THEN 
DO 201 N ■ 1 , IREC 
XI ( N ) - X< I,N) 

YI(N) - H( I , N) 

CONTINUE 

ELSE 

DO 202 N « 1, IREC 
NM - IREC + 1 - N 
XI ( N ) - X( I , NM) 

YI(N) - H { I , NM) 

CONTINUE 
END IF 
L - 1 

DO 203 N * 2, IREC 

IF ( XI (N). GT. XI ( 1 ) ) THEN 
L » L + 1 
XI (L) ■ XI (N) 

YI ( L ) - YI(N) 

END IF 
CONTINUE 
IRECA * IREC 
IREC - L 

WRITE ( 6 , * ) ' ' 

WRITE ( 6,* ) ' ' 

WRITE ( 6,* ) ' 9 , 1 SPLINE FITTING TO PASS : ' , NREV< I ) , ' f 

IREC , * DATA POINTS' 

IF ( IRECA. NE. IREC ) THEN 

NJ - IRECA - IREC 

WRITE ( OUP, * ) ' ' , N J , ' DATA POINTS THAT ARE NOT IN ', 

'ASCENDING OR DECENDING ORDER HAVE BEEN REMOVED.' 

END IF 

CALL ICSCCU ( XI, YI, IREC, CIMSL, IREC-1, IER ) 

IF ( IER. NE. 0 ) THEN 

WRITE < OUP, 609 ) (J,XI(J) , YI ( J ) ,J*1, IREC ) 

END IF 
L » 0 

DO 204 J - 1, NPASS 

IF ( ICX(I,J). NE. 0 ) THEN 
L = L + 1 
U(L) » EX( I , J) 

WRITE ( OUP, * ) ' PASS: ' ,NREV( J) , ' ELON = 

EX( I , J) , ' U = ' ,U(L) 

END IF 
CONTINUE 

WRITE ( OUP, * ) ' NUMBER OF CROSSOVERS = ',L 

IF ( L. GT. 0 ) THEN 

CALL ICSEVU ( XI, YI, IREC, CIMSL, IREC-1, U, S, L, IER ) 
IF ( IER. NE. 0 ) THEN 

WRITE ( 6 , * ) ' ' , ’ ELON ARRAY ' 

WRITE ( 6,* ) (XI(MM) ,MM=1, IREC) 

WRITE ( 6,* ) ' HSS ARRAY' 

WRITE ( 6,* ) ( YI( MM) , MM= 1 , IREC) 

END IF 
L = 0 

DO 205 J = 1, NPASS 

IF ( ICX(I,J). NE. 0 ) THEN 
L = L + 1 
XVR( I , J ) = S(L) 

WRITE ( OUP, XPASS : ' , NREV ( J ) , ' ELON 

U(L) , ' HSS = ' ,S<L) 

END IF 
CONTINUE 
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END IF 
200 CONTINUE 


NPASSl - NPASS - 1 
WRITE ( OUP, 600 ) 
WRITE ( OUP, * ) ' 

WRITE ( OUP, 608 ) 


' , ' SUMMARY OF CROSSOVERS DIFFERENCES 
^BEFORE ADJUSTING FOR BIASES : ' 


CALL XOVD ( NPASS, IDS, ICX, ISUM, NREV, MAXOVR, XVR AX BX 
‘ EX, AXTBX, XDMAX, KC, XRMS, ' ' BX 

4 W, AXT, SMSS, SMGG, SMSG ) 


DO 300 I m 1 
DO 301 J * 
AX( I , J- 1 ) 
301 CONTINUE 
300 CONTINUE 


KC 

IXM+1, NPASS 
- AX(I,J) 


DO 400 I - 
DO 401 J 
AXT( J, I) 
401 CONTINUE 
400 CONTINUE 


1, KC 

- 1, NPASSl 
- AX( I, J) * 


W<I> 


^CALL VMULFF ( AXT^AX, NPASSl, KC, NPASSl, MPASS, 
CALL VCVTFS ( SCR, NPASS 1 \ MPASS, AXTAX ) 


MAXOVR, SCR, 


CALL VMULFF ( AXT, BX, NPASSl, KC, 
4 MPASS, IER ) 


MPASS, MAXOVR, AXTBX, 


DO 700 I a 
IDIAG( I) 

700 CONTINUE 

DO 701 I - 
DAXTAX ( I ) 

701 CONTINUE 

DO 702 I - 
DAXTBX ( I ) 

702 CONTINUE 


1, NPASSl 

1 ( I * (1+1) ) / 2 

1, IDIAG( NPASSl) 

■ DBLE ( AXTAX ( I ) ) 

1, NPASSl 
- DBLE ( AXTBX ( I ) ) 


CALL SOLVE ( DAXTAX, DAXTBX, ID I AG, 


NPASSl, 


DO 800 I - 
K - NPASS 
DAXTBX ( K ) 
800 CONTINUE 

DAXTBX ( IXM) 


IXM, NPASSl 
+ IXM - I 
» DAXTBX (K-l) 


. TRUE . , 


. TRUE . ) 


DO 801 I 
AXTBX ( I) 
801 CONTINUE 


1 , NPASS 
SNGL< DAXTBX < I) ) 


XRB = 0. 

DO 802 I a 
XRB * XRB 
802 CONTINUE 

XRB = XRB / 


1, NIREC ( IXM ) 

+ H( IXM, I ) - HG ( IXM , I ) 

FLOAT ( NIREC ( IXM) ) 


DO 803 I 
AXTBX ( I ) 
803 CONTINUE 


1 , NPASS 
AXTBX ( I ) + XRB 


WRITE ( OUP, 600 ) 

WRITE ( OUP, * ) ' 

WRITE ( OUP, 608 ) 


' SUMMARY OF CROSSOVERS DIFFERENCES 
AFTER ADJUSTING FOR BIASES:' 

XDMAX 


. XOVD ( NPASS, IDS, ICX, ISUM, NREV, MAXOVR, XVR AX BX 

1 EX, AXTBX, 1000.0, KC, XRMsI ' ' X ' 



non 


* W, AXT, SMSS, SMGG, SMSG ) 

C 

WRITE ( OUP, 601 ) NREV(IXM), XRB 
DO 900 1*1, NPASS 

IF ( ISUM(I). EQ. 0 ) THEN 
AXTBX(I) *0. 

WRITE ( OOP, 603 ) I, NREV(I), AXTBX(I) 

C ELSE IF < I. EQ. IXM ) THEN 

C WRITE ( OUP, 604 ) IXM, NREV( IXM) 

ELSE 

WRITE ( OUP, 602 ) I, NREV { I ) , AXTBX(I) 

END IF 
900 CONTINUE 

DO 1100 1*1, NPASS 

ICB(I) * IFIX(AXTBX( I) *1000. ) 

1100 CONTINUE 

REWIND INP 
NREC * 0 

4000 READ ( INP, END - 9000 ) II , 12, 13, 14, IS, 16, R1,R2 , 17 , 18, R3, 

‘ 19, I 10, 111, R4, 112, 113, I 14, 115, I 16, 117, 118, I 19, 

6 120,121,122,123,124,125,126,127,128 


IF ( 114. EQ. 12793 ) 
113 ■ -1000000 
GO TO 888 
END IF 


THEN 


NREC = 
IB IAS 


NREC + 1 

- 0 


DO 1101 I » 1, NPASS 

IF ( 114. EQ. NREV ( I ) ) THEN 

I8IAS * I 
GO TO 5000 
END IF 
1101 CONTINUE 
C 

5000 CONTINUE 
C 

IF ( IBIAS . EQ. 0 ) THEN 
113 * -1000000 
ELSE 

113 - 113 - ICB(IBIAS) 

END IF 

888 CONTINUE 

WRITE ( OUG ) I1,I2,I3,I4,I5,I6,R1,R2,I7,I8,R3, 

‘ 19, 110, 111, R4, 112, 113, 114, 115, 116, 117, 118, 119, 

* ~~ 120, 121, 122, I2j, 124, 125, 126, 127, 128 

GO TO 4000 

9000 WRITE ( OUP, 600 ) 

WRITE < OUP, * ) ' NREC, ' RECORDS ADJUSTED FOR ORBIT BIAS’ 

DETERMINE THE ERROR COVARIANCE MATRIX 

CALL SYMINV ( DAXTAX , DAXTBX , ID I AG, NPASS 1, ERCOVM ) 

MM = ( NPASS 1* ( NPASS 1+ 1 ) ) / 2 

WRITE ( OUC, 610 ) { ERCOVM ( I ) , I = 1 , MM ) 

WRITE { OUP, * ) ' ' 

WRITE { OUP, * ) ' 


'THE COVARIANCE MATRIX OF BIAS ADJUSTMENT' 
' HAS BEEN WRITTEN TO FILE BLACK. PCOV.' 


501 FORMAT ( 2X , 1 3 , 8F9 . 3 ) 

600 FORMAT ( 1H1 ) 

601 FORMAT( 1H1,T20, 'BIAS FOR SATELLITE PASSES OVER THE 
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in J.U) ) 

607 FORMAT((lH , 10 ( 1PE12 . 4 ) ) ) 

608 FORMAT(lH , /,T10,'THE CUTOFF WINDOW FOR CROSSOVER DIFFERENCES' 

„ 6 ' IS ' , F7 . 2 , ' METERS ' / ) ' 

609 FORMAT ( ( 1H , T1 0 , 16 , 2F12 . 4 ) ) ' 

610 FORMAT ( (4D20.13) ) 


STOP 

END 


SUBROUTINE XOVD ( 


NPASS , IDS, ICX, I SUM, NREV, 
AX, BX, EX, BIAS, XDMAX, KC, 
W, AXT , SMSS, SMGG, SMSG ) 


MAXOVR, 
XRMS , 


XVR, 


REAL* 4 XVR ( NPASS , NPASS ) , EX ( NPASS , NPASS ) 

REAL 4 BIAS ( NPASS ) , AX ( MAXOVR, NPASS) , BX(MAXOVR) 

REAL*4 W( MAXOVR), AXT ( NPASS , MAXOVR ) ' 

£?2 <NP ?f S ' M £e SS) ' IDS (NPASS), ISUM( NPASS) , 
INTEGER* 4 K75, K78, K7578, KC, MAXOVR, OUP 
DATA OUP / 6 / 

XRMS * 0. 

KC - 0 
K75 - 0 
K78 - 0 
K7578 - 0 

DO 100 I » 1, NPASS 
DO 101 J - I, NPASS 

IT ( ICX(I,J). NE. 0 ) THEN 
XVR( I , J) - XVR ( I , J ) - BIAS ( I ) 

XVR ( J , I ) * XVR ( J , I ) - BIAS ( J ) 

BXT - XVR( I , J) - XVR( J, I ) 
rry it t * irr* i * 


NREV ( NPASS j 


IF ( ICX(I,J) 
ICX(I,J) - 0 
ICX(J,I) - 0 
WRITE ( OUP , 


EQ. -1 ) THEN 


600 ) 


EX( I, J) , NREV ( I ) , XVR( I, J) , NREV ( J ) , 
XVR(J,I),BXT 
LE . XDMAX ) THEN 


ELSE IF ( ABS(BXT) 

KC » KC + 1 
AX ( KC , I ) » 1.0 
AX(KC,J) - -1.0 
BX(KC) = BXT 
XRMS ■ XRMS + BXT* BXT 

WRITE ( OUP, 601 ) EX( I, J) , NREV( I ) , XVR( I, J) NREV( J) 

XVR ( J , I ) , BX ( KC ) 

IT ( IDS < I ) . EQ. IDS ( J ) ) THEN 

IF ( IDS (I) . EQ. 7502701 ) THEN 
1 


K75 = K75 
W ( KC ) * SMGG 
END IF 

IF ( IDS ( I ) . EQ. 


K78 
W ( KC ) 
END IF 
ELSE 

K7578 = 
W(KC) = 
END IF 
ELSE 

ICX< I, J) 
ICX ( J , I ) 


K78 + 1 
= SMSS 


7806401 ) THEN 


K7578 

SMSG 


= 0 
= 0 


+ 1 
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u u u u 


XSUM(I) 

» ISUM(I) 

ISUM( J ) 

- ISUM(J) 

WRITE ( 

OUP, 602 

& 


END IF 
END IF 


101 CONTINUE 

100 CONTINUE 



1 

1 

EX ( I , J ) , NREV ( I ) 
XVR ( J , I ) , BXT 


XVR( I, J) , NREV ( J ) , 


C 

C 

C 


C 


XRMS - SQRT ( XRMS / FLOAT (KC) ) 


WRITE ( OUP, 603 ) XRMS, K75, K78, K7578, KC 


600 FORMAT ( 1H 


601 


602 


603 


FORMAT ( 1H 


FORMAT ( 1H 


FORMAT ( 1H0 


, T15 , ' X § E. LONG - ' , F6 . 2 , 8X, 

,15, ' ( ' ,F6.2, ' ) - #' ,15, ' ( ',F6.2, ' ) ' , 

' > ' , F8 . 2 , 5X, ' ** EDITTED OUT (A PRIORI) 

, T15 , ' X § E. LONG - ',F6.2,8X, 

' # ' * 15, ' ( ' , F6 . 2 , ' ) - # ' , 15, ' ( ' , F6 . 2 , ' ) ' , 

' > ' , F8 . 2 ) 

,T15,'X 0 E. LONG - ' , F6 . 2 , 8X, 

'#' ,15, ' ( ' ,F6.2, ' ) - #' ,15, ' ( ' ,F6.2, ' ) ' , 

' > ' , F8 . 2 , 5X, ' ** DISCARDED (IN SITU) **' 

,/,T20,'RMS (M) - ' , F10 . 4, // , 

T24 , ' # OF GE03-GE03 CROSSOVERS = ',15,//, 

T24 , ' # OF SEASAT-SEASAT CROSSOVERS - ',15,//, 

T24 , ' I OF GE03-SEASAT CROSSOVERS - ',15,//, 

T2 4, 'TOTAL # OF CONSIDERED CROSSOVERS =■ ',15) 


** > 


) 


C 


c 


RETURN 

END 


SUBROUTINE CURVE ( NPASS, X, Y, Z, ZG, YR, IREC, IDP, IDSAT, 
‘ A, B, C, F, FT, G, 

4 SX, SY, SZ, SZG, DDZ , D, SIGE ) 


REAL* 4 
REAL* 4 
REAL* 4 
REAL* 8 
REAL* 8 
REAL* 8 
REAL* 8 
INTEGER* 4 
INTEGER* 4 
INTEGER* 4 


X ( IREC ) , Y(IREC), Z(IREC), ZG(IREC), YR(IREC) 
SX(IREC), SY(IREC), SZ(IREC) 

SZG (IREC), DDZ (IREC), D(IREC) 

F ( IREC , 3 ) , FT ( 3 , IREC ) 

E<3,3), El ( 3 , 3 ) , AE ( 6 ) 

G(IREC), H<3) 

A, B, C, DET, YD 

IDP, IDSAT, IREC, NPASS 

IDIAG<3) / 1, 3, 6 / 

OUP / 6 / 


EVICT THE BAD DATA POINTS BY COMPARING THE ABSOLUTE DIFFERENCE 
BETWEEN ADJACENT DATA POINTS 


NB ■ 0 

IREC 1 = IREC - 1 

C 

DO 111 1=1, IREC 

SX( I) = -1.0 

111 CONTINUE 
C 

DO 112 1=1, IREC 1 

DDZ ( I ) = ABS ( Z(I+1) - Z < I ) ) 

112 CONTINUE 
C 

DO 113 J = 1, IREC 1 

IF ( DDZ ( J ) .LT.SIGE) THEN 
K = J 
GO TO 200 
END IF 

113 CONTINUE 
C 

200 CONTINUE 

DO 222 J = 1, IREC 

IF(DDZ( IREC-J) .LT.SIGE) THEN 
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c 


c 


IRECG - IREC-J 
GO TO 300 
END IF 
222 CONTINUE 


300 KB - 0 

DO 333 I - 1, IRECG- 1 
ir(DDZ( I) . GT.SIGE) 
333 CONTINUE 


KB 


KB + 1 


C 


C 


LG - 1 

SX(LG) * X ( K ) 
SY(LG) - Y(K) 
SZ(LG) - Z(K) 
SZG(LG) * ZG(K) 


400 LG * LG + 1 
K - K + 1 

IF ( DDZ ( K— 1 ) .GT.SIGE) 
K - K + 1 
DO 444 I 


THEN 


1, KB- 1 

IF ( DDZ ( K-l ) .GT.SIGE) 
444 CONTINUE 
END IF 

SX(LG) - X(K) 

SY(LG) * Y(K) 

SZ(LG) - Z(K) 

SZG(LG) - ZG( K ) 

IF( K. LE. IRECG) GO TO 400 


K + 1 


DO 555 1*1, LG 

F( I, 1) ■ SX( I) *SX( I) 
F( 1 , 2 ) - SX< I ) 

F(I,3) - 1.0D0 
G(I) * SY ( I ) 

555 CONTINUE 


DO 666 I * 
DO 667 J 
FT ( J , I ) 

667 CONTINUE 
666 CONTINUE 


1, LG 
■1,3 
■ F(I,J) 


C 


DO 777 1*1 

DO 770 J - 
E( I, J) * 
DO 779 K 
E( I, J) i 

779 CONTINUE 

778 CONTINUE 
777 CONTINUE 


, 3 
1, 3 
0 . DO 
» 1, LG 
* E ( I , J ) 


+ 


FT( I , K) *F ( K, J) 


C 


C 


DO 888 I » 1, 3 

H(I) = 0 . DO 
DO 889 K * 1, 

QOO H<1) = H < Z > 
889 CONTINUE 

888 CONTINUE 


LG 

FT ( I , K ) *G 


(K> 


AE(1) - E ( 1 , 1 ) 
AE<2) * E( 1,2) 
AE ( 3 ) = E ( 2 , 2 ) 
AE ( 4 ) = E ( 1,3) 
AE ( 5 ) = E( 2 , 3 ) 
AE ( 6 ) - E ( 3 , 3 ) 


C 


CALL SOLVE ( AE , H, 


IDIAG, 3, .TRUE., . 


A = H ( 1 ) 
B = H( 2 ) 


TRUE. 
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C - H<3) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


DET - E<1,1) * ( E(2,2)*E(3,3) 
4 E( 1,2) * ( E(1,3)*E(2,3) 
6 E( 1,3) * ( E ( 1,2) *E (2,3) 


E(2,3)*E(2,3) ) +■ 
E ( 1 / 2 ) *E( 3,3) ) + 
E(1,3)*B<2,2> ) 


EI(1,1) - ( E(2,2)*E(3,3) 
EI<2,2) - ( E(1,1)*E(3,3) 
El ( 3 , 3 ) - ( E ( 1 , 1 ) *E( 2,2) 
El < 2 , 3 ) - ( E( 1 , 2 ) *E( 1,3) 
EI(1,3) - ( E(1,2)*E<2,3) 
EI(1,2) - ( E ( 1,3) *E (2,3) 


- E( 2 , 3 ) *E( 2,3) ) / DET 

- E(1,3)*E(1,3) ) / DET 

- E( 1,2)*E< 1,2) ) / DET 

- E ( 1 , 1 ) *E (2,3) ) / DET 

- E( 1,3)*E(2,2) ) / DET 

- E( 1,2)*E(3,3) ) / DET 


A - ( El ( 1 , 1 ) *H ( 1 ) + El ( 1 , 2 ) *H( 2 ) + El ( 1 , 3 ) *H( 3 ) > 

B - ( El ( 1 , 2 ) *H( 1 ) + El ( 2 , 2 ) *H( 2 ) ♦ EI(2,3)*H(3) ) 

C * ( EI( 1, 3 ) *H( 1 ) + EI( 2 , 3 ) *H( 2 ) + EI(3,3)*H(3) ) 

RMS » 0. 

DO 999 I ■ 1, LG 

YD - A*SX(I)*SX(I) + B*SX ( I ) + C 
YR( I ) - SNGL( YD) 

D( I ) - YR( I ) - SY ( I ) 

RMS * RMS + D( I ) *D( I ) 

999 CONTINUE 

RMS - SQRT ( RMS / FLOAT ( LG ) ) 


C 


WRITE ( OUP, 601 ) 

4 

IT ( RMS. GT. 0.005 

WRITE < 6, * ) 

WRITE ( 6, 602 ) 

WRITE ( 6, * ) 

END IF 


NPASS , IDT, IDSAT, A, B f C, 

SX(1), SY(1), SX(LG), SY (LG) , RMS, 
) THEN 

t 

(M,SX(M) ,SY<M) ,YR(M) ,D<M) ,M-1,LG) 


IF { LG. NE. IREC ) THEN > 

IREC » LG 

DO 1111 I-l, IREC 
X(I) - SX ( I ) 

Y(I) - SY ( I ) 

2(1) » SZ<I) 

ZG(I) » SZG(I) 

1111 CONTINUE 
END IF 

601 FORMAT ( 1H0, // , IX, 13 , 17 , 19 , 3E14 . 5 , 5F10 . 4 f 16 ) 

602 FORMAT ( ( 1H , T5 , I 4 , 3X , 3F 1 0 . 4 , 2X , F 1 0 . 4 ) ) 


RETURN 

END 


LG 


C 


C 


SUBROUTINE XOVER ( NPASS, NREV, IDSAT, A, B, C, SLONl, SLATl, 

& SLON2, SLAT2 , ICX, EX, GY, ISUM, LN, LY, IXM ) 


INTEGER 

INTEGER 

INTEGER 

REAL *4 

REAL *4 

REAL *4 

REAL* 8 

REAL* 8 

REAL *8 

INTEGER 

INTEGER 

INTEGER 

DATA 

DATA 

DATA 


NREV ( NPASS ) , IDSAT ( NPASS ) 

ICX( NPASS, NPASS) , ISUM( NPASS), LN< NPASS), LY( NPASS) 
INDEX ( 2 ) ' 

SLATl { NPASS ) , SLAT2 ( NPASS ) , SLON 1 ( NPASS ) , SLON2 ( NPASS ) 
EX( NPASS, NPASS) , GY < NPASS , NPASS ) 

P , Q 

A( NPASS), B( NPASS), C( NPASS) 

T ( 2 ) 

DA, DB, DC, RT, SRT 

ITOTAL, IXM, MAX IX , MX, MXS , MX75, MX78, NPASS 
IJK, IP, JP 
JOUT, OUP 

JOUT / 4 /, OUP / 6 / 

MAXIX / 0 /, MX / 0 /, MX7 5 / 0 /, MX78 / 0 / 

ITOTAL / 0 / 


DO 123 1=1, NPASS 
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u u 


ISUM(I) » 0 
DO 234 J - 1, NPASS 
ICX(I,J) « 0 
EX(I,J) - 0. 

GY ( I , J ) - 0. 

234 CONTINUE 
123 CONTINUE 
C 

WRITE ( OUP, 600 ) 

C 

DO 222 I - 1, NPASS- 1 

c 

WRITE ( OUP, 601 ) NREV(I), IDSAT(I), SLONl(I), SLATl(I) 
c 6 SLON2 ( I ) , SLAT 2 ( I ) ' ' ' 

DO 333 J ■ 1+1, NPASS 
C 

INDEX ( 1 ) * 0 
INDEX( 2 ) » 0 
DB - B ( J) - B ( I ) 

DA ■ A( I ) - A(J) 

DC - C(I) •- C)J) 

RT ■ DB*DB - 4 . 0*DA*DC 

C 

IF ( RT. GE. 0.0 ) THEN 
IF ( RT. EQ. 0.0 ) THEN 
SRT - 0.0 
ELSE 

SRT * DSQRT(RT) 

END IF 

T( 1 ) ■ ( DB+SRT) / ( 2 . 0*DA) 

T( 2 ) - ( DB-SRT) / ( 2 . 0*DA) 

WRITE ( 6 , * ) T, RT, SRT 


4 


4 


4 


4 


135 

C 


4 


246 


DO 135 L * 1, 2 
IF ( SLON2 ( I ) . 
IF ( T(L). GE. 
INDEX(L) - 

ELSE 

IF ( T(L). LE. 
INDEX ( L) - 
END IF 

IF ( SLON2 ( J ) . 
IF( T(L). GE. 
INDEX ( L) » 

ELSE 

IF( T(L) . LE. 
INDEX ( L ) * 
END IF 
CONTINUE 


GT. SLON 1(1) ) THEN 

SLON 1(1). AND. T(L) . LE. 
> INDEX ( L ) + 1 

SLONl(I). AND. T(L). GE. 
1 INDEX) L) + 1 

GT. SLON 1 ( J ) ) THEN 

SLONl(J). AND. T(L) . LE. 
INDEX) L) +■ 1 

SLON 1 ( J ) . AND. T) L) . GE. 
INDEX) I,) + 1 


SLON2 ( I ) 


SLON2 < I ) 


SLON 2 ( J ) 


SLON 2 ( J) 


) 

) 


) 

) 


IF ( INDEX) L) 

MX ■ MX + 1 
P - SNGL ( T ( L ) ) 

Q * SNGL ( A) I ) *P*P 


EQ. 2 ) THEN 


EX) I, J) 
GY(I,J) 
EX) J, I) 
GY ( J , I ) 
WRITE ( 


= P 
■ 0 
= p 
= Q 
OUP, 


+ B( I ) *P + C 


602 ) 


NREV(J), iDSfl 
SLON 1 ( J ) , SLAT 
ICX(I,J) = ICX ( I , J ) + 1 
ICX(J,I) = ICX) I, J) 

IF ( IDSAT(I). EQ. IDSAT(J) 

IF ( IDSAT(I). EQ. 7502701 


) 


IF ( 
END IF 
END IF 
CONTINUE 


IDSAT(I). EQ. 7806401 


) 


. P, Q, 

, SLON 2 ( J ) , SLAT2 ( J ) 


= MX 7 5 + 1 
= MX7 8 + 1 
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non 


c 

END IF 
C 

333 CONTINUE 
C 

222 CONTINUE 


WRITE 

< 

OUP, 

* ) 

t 

r 





WRITE 

( 

OUP, 

* ) 

t 

t 





WRITE 

( 

OUP, 

* ) 

t 

TOTAL 

* 

or 

POSSIBLE CROSSOVERS - ' , 

MX 

WRITE 

( 

OUP, 

* ) 

* 

f 





WRITE 

( 

OUP, 

* ) 

t 

TOTAL 

» 

or 

GE03-GE03 CROSSOVERS * ' , 

MX75 

WRITE 

( 

OUP, 

* ) 

r 

* 




WRITE 

( 

OUP, 

* ) 

f 

TOTAL 

# 

or 

SEASAT-SEASAT CROSSOVERS 

- ', MX7 8 

MXS - 

MX - MX 7 5 

- MX7 8 





WRITE 

( 

OUP, 

* ) 

/ 

f 





WRITE 

( 

OUP, 

* ) 

* 

TOTAL 

* 

or 

GEOS 3 -SEAS AT CROSSOVERS * 

' , MXS 


EDIT OUT BAD CROSSOVERS BASED ON A PRIORI KNOWLEDGE 


DO 555 IR * I, 500 

READ ( JOUT, 501, END = 556 ) JC, IP, JP 
IF ( JC. EQ. 0. AND. ICX(IP,JP). NE. 0 ) THEN 
ICX(IP,JP) * -1 
ICX(JP,IP) - -1 
END IF 

555 CONTINUE 
C 

556 CONTINUE 
C 

DO 345 1*1, NPASS 

DO 456 J - 1, NPASS 
IJK • ICX(I,J) 

IF ( IJK. LT. 0 ) IJK - 0 
ISUM(I) ■ ISUM(I) + IJK 
456 CONTINUE 

ITOTAL - I TOTAL ♦ ISUM(I) 

345 CONTINUE 

ITOTAL * ITOTAL/ 2 
C 

DO 567 I - 1, NPASS 

MAX IX * MAXO ( MAXIX, ISUM(I) ) 

567 CONTINUE 


C 


C 


C 


C 


DO 678 1*1, NPASS 

IF ( ISUM(I). EQ. MAXIX ) IXM * I 
678 CONTINUE 

WRITE ( OUP, 603 ) 

WRITE ( OUP, 604 ) (I, NREV(I), 1=1, NPASS ) 

WRITE ( OUP, 605 ) NPASS 

WRITE ( OUP, 606 ) ( NREV(I) , ( ICX( I , J) , J*1 , 94 ) , ISUM< I ) , 

*■ 1*1, NPASS ) 

WRITE ( OUP, 607 ) NREV(IXM), MAXIX, ITOTAL 

MM = 0 
NN * 0 

DO 789 1=1, NPASS 

IF ( ICX( IXM, I ) . EQ. 1 ) THEN 
MM = MM + 1 
LY(MM) = I 
ELSE 

NN = NN + 1 
LN ( NN ) = I 
END IF 
789 CONTINUE 

WRITE ( OUP, 608 ) NREV(IXM), MM, NN, MM+NN 
WRITE ( OUP, 609 ) ( NREV ( LY ( I ) ) , I = 1 , MM ) 
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oononooonooonnooonon 


876 

888 


DO 987 J « 1, NN 

WRITE ( OUP, 610 ) NREV ( LN ( J ) ) 

DO 876 I - 1, NPASS 

ir ( ICX( I , LN ( J) ) . EQ. 1 ) THEN 
WRITE ( OUP, 611 ) NREV(I) 

GO TO 888 
END IF 
CONTINUE 
CONTINUE 


987 CONTINUE 
501 FORMAT ( 11,2110) 

603^FORMAT< 1H1,T20, 'ROW-COLUMN ELEMENTS OF PASS CORRELATION MATRIX:', 
604 FORMAT ( ( 1H , 10< 3X, 14, 16) ) ) 

6°5/ORMAT(lHl,T20, 'PARTIAL (94) PASS CORRELATION MATRIX FOR', 

607 FORMAT( lHO,//?Tiof ; PASS 4 | 1; ?i6“ , Ls THE MAXIUM OF ' 15 
® 08 FORMAT ( lH??f f ' I H =, ^£iS X CROSSOVERS ' IS : ' , IS ) 

«o. 4 «»uT„i H oTi!i;r DIRECT cross4as «s Arif?:/'/, 

FORMAT ( TlO^REFERENCE * PASS v£ T ° ™ E ' > 

RETURN 

END 


SUBROUTINE SOLVE (A, B, IDIAG, NEQ, FACT,. BACK) 


S3£d"i?2 U n T FAC? I ™uE aCt ^i Zat t° n ° f the s y"®»tric matrix 

in a, ir FACT - TRUE; and solve A * X - B if BACK - TRUE. 

A oArf a «?“*.K h ® co ??? c ^ ed_column form of the upper trianqular 

Addresses of the diagonal terms in A. 

Number of equations 

Tf n irl * th ? n factor A ? otherwise do not factor A 

do b “ l '»“ b «'itut.; otherwise 


IDIAG 

NEQ 

FACT 

BACK 


IMPLICIT REAL* 8 (A-H.O-Z) 

LOGICAL FACT, BACK 

DIMENSION A( 1 ) , B ( 1 ) , IDIAG ( 1 ) 


Factor A, reduce B 
JR - 0 

DO 400 J » 1 , NEQ 
JD = IDIAG ( J ) 

JH = JD - JR 

IS =* J - jh + 2 

IF (JH . LT. 2) GOTO 390 
IF (FACT) THEN 
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IF (JH .GT. 2) THEN 

C. . Reduce column J rows IS to J— Is do not divide by row diagonal 

K * JR + 2 
ID « IDIAG ( IS - 1) 

C 

DO 100 I * IS, J-l 
IR * ID 
ID » IDIAG ( I ) 

IH * MIN ( ID-IR-1, I-IS+1) 

IF (IH.GT.O) A(K) - A(K) - DOT( A( K-IH ) , A( ID-IH ) , IH ) 

K - K + 1 

100 CONTINUE 

ENDIF 
C 

Divide by row diagonal , and reduce diagonal term in column J 
C 

IR * JR + 1 
K - J - JD 

DO 200 I - IR, JD- 1 
ID * IDIAG ( K+I ) 

IF ( A{ ID ) • EQ« 0.0) GOTO 200 
D - -A(I) 

A<I) - A( X) /A ( ID ) 

A( JD ) * A< JD ) + D*A( I ) 

200 CONTINUE 

C 

ENDIF 

C 

C • « Reduce RHS 

IF (BACK) B ( J ) =■ B ( J ) - DOT{ A( JR+1 ) , B( IS— 1 ) , JH- 1 ) 

C 

390 JR ■ JD 
400 CONTINUE 
C 
C 

IF ( . NOT . BACK ) RETURN 
C 
C 

C Divide by diagonal pivots 


DO 700 I » 1 , NEQ 
ID * IDIAG ( I ) 

IF ( A( ID ) . NE .0.0) B ( I ) 
700 CONTINUE 
C 

C Backsubstitute 
C 

J =* NEQ 

JD =* IDIAG ( J ) 

C 

801 D * -B ( J ) 

J - J - 1 

IF (J.LE.0) RETURN 
C 


B < I ) / A( ID ) 


JR = IDIAG ( J ) 

IF (JD-JR.GT.l) THEN 
IS =* J - JD + JR + 2 
K > JR - IS M 
DO 810 I = IS, J 
810 B ( I ) =* B ( I ) + A(I+K)*D 
ENDIF 
C 

JD = JR 
GOTO 801 
C 

END 

C 

FUNCTION DOT <A,B,N) 

C 

C Compute the dot product of the two N-vectors A and B. 
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c 

c 


INTEGER* 4 
REAL *8 

DOT * 0.0 


N, I 

DOT, A(1),B(1) 


C 

C 


c 

c 


DO 100 I - 1, N 

DOT - DOT + A( I ) *B ( I ) 
100 CONTINUE 

RETURN 

END 

SUBROUTINE CLR { A , N A ) 

INTEGER* 4 NA, I 
REAL*4 A(l) 

DO 100 I - 1, NA 
A(I)»0.0 
100 CONTINUE 

RETURN 

END 

SUBROUTINE DCLR ( A , N A ) 


INTEGER* 4 
REAL* 8 


NA, I 
A(l) 


.FALSE., .TRUE. ) 


DO 100 I - - 1, NA 

A( I)-0. 0D0 
100 CONTINUE 
C 

RETURN 
END 

c SUBROUTINE SYMINV ( A, B, ID1AG, NEQ, AINV ) 

IMPLICIT REAL* 8 (A-H,0-3) 

c DIMENSION A(l), B<1), AINV(l), IDIAG(l) 

K - 0 

DO 100 1-1, NEQ 

CALL DCLR ( B, NEQ ) 

B ( I ) - 1.0 

CALL SOLVE ( A, B, IDIAG, NEQ, 

DO 110 J » 1, I 
K a K + 1 
AINV(K) - B(J) 

110 CONTINUE 
100 CONTINUE 
C 

RETURN 
END 

//* 

// EXEC LINKGOV, REGION=5000K 
//SYSLIB DD DSN=SYS2 . IMSLS , DISP=SHR 
DD DSNa S YS2.IMSLD,DISP=SHR 

//★ *** + ***★ **** * + ** **** **** 

//* 

//* 

//* 

//* **** *** + ** + ★ **** **** *** 

//* 

/ /22’ DD DSN=zm AYA. ALTIM. DATA( BLKRAPP ) DISP-SHR 

DD d sn-zmaya. black. data; disp=shr 
//GO. FT04F00 1 DD DSN=ZMAYA. BLKPOUT . DATA DISP=SHR 
//GO.FT06F001 DD SYSOUT=* ' DISP “ SHR 

//GO. FT08F001 DD DSN=ZMAYA. BLACKA. DATA, DISP=SHR 


* * * * 


* ★ ★ ★ ★ * * * * 


BLACK SEA DATA 


*★★★ *♦++ ★ ★ * i 
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//GO.FT09F001 DD DSN-ZMAYA. BLACK. PCOV, DISP-SHR 
// EXEC NOTIFYTS 
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on o 


A listing of program PLTAEP 


/ * JOBPARM A LINES-6 0 109 ' 360 ' 75) ' AYAU .TIME- (12,00), CLASS-O , MSGCLASS-X 

//TEMPLATE EXEC G38PLOT 
/ /FORT . SYS IN DD * 

C 

C rORMAT OP GEOS-3/SEASAT ALTIMETER DATA 


DESCRIPTION 
SATELITE ID 

MEASUREMENT TYPE ( 42- OVER LAND , - 
TIME SYSTEM ( NM ) 

STATION NUMBER 

PREPROCESSING INDICATORS 

MODIFIED JULIAN DATE OF OBSERVATION 


OVER WATER ) 


FRACTION OF DAY PAST MIDNIGHT (GMT) 

ALTIMETER OBSERVATION ( METERS ) 

SATELLITE GEODETIC LATITUDE ( IE-6 DEGREES) 

SATELLITE EAST LONGITUDE ( IE— 6 DEGREES ) 

MEASUREMENT STANDARD DEVIATION (METERS) 

NET INSTRUMENT CORRECTION (MM) 

METEOROLOGICAL DATA WORD (GEODYN VOL 3) 

NET MEDIA CORRECTIONS (MM) 

GEOID HEIGHT ABOVE REFERENCE ELLIPSOID (METERS) 
NET OCEAN DYNAMIC CORRECTIONS (MM) 

INDICATED SURFACE ELEVATION (MM) 

S/C REVOLUTION NUMBER 

MEAN SEA SURFACE ELEVATION (MARSH/MARTIN '81 (MM)) 
DOD REFERENCE RADIAL ORBIT. DIFFERENCE (MM) 

ag £ /3 (CM) 

AGC ( OB ) 

WIND SPEED (' CM/SEC) 

SURFACE ELEVATION PREPROCESSING WORD 
DRY TROPOSPHERIC CORRECTION (MM) 

FNOC WET TROPOSPHERIC CORRECTION (MM) 

SMMR WET TROPOSPHERIC CORRECTION (MM) 

IONOSPHERIC CORRECTION (MM) 

BAROTROPIC DYNAMIC SEA SURFACE CORRECTION (MM) 
SOLID EARTH TIDE ( MM) 

SCHWIDERSKI OCEAN TIDE (MM) 

PARKE OCEAN TIDE MM) 


PARAMETER ( MAXPLT 
PARAMETER ( NMAX - 
PARAMETER ( NINT - 
PARAMETER ( ETOP - 
PARAMETER ( EBOT - 
PARAMETER ( SIGE - 
PARAMETER ( SIGT - 


- 9999 ) 
2500 ) 
150 ) 
100 . ) 
0 . ) 
2.0 ) 
8.0E-5 ) 


PARAMETER { IC * NMAX - 1 ) 


INTEGERS 
INTEGER* 4 
REAL* 8 
REAL *4 
INTEGER* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
INTEGER 
INTEGER 
INTEGER 


Jq 7 ;! 1 ^ 1 ^'^ 0 '??!' 122 ' 123 ' 124 . 125,126,127,128 

R 1 , R2 
R3 , R4 
IDSAT 

GLAT, ELON, HSS 

X( NMAX) , Y ( NMAX ) , YM( NMAX ) , DF ( NMAX) , C{ NMAX-1 , 3 ) 

F ( NMAX ) , GX ( NMAX ) , GY ( NMAX ) , DGY ( NMAX ) , DD Y ( NMAX ) 

DDT (NMAX ) , XA( NMAX ) , YA( NMAX ) , DYA(NMAX) 

WK ( 7 *NMAX+ 14), XF(NINT+1), YF ( NINT+ 1 ) 

INP, IREC, I PLOT , NREC, NREV, OUP 
LIST ( 6 ) 

IBC(NMAX) 


DATA INP / 3 /, ING / 4 /, OUP / 6 / 

DATA LIST / 8, 110, 75, 0, 7, 10 / 

DATA IREC / 0 /, NREC / 0 /, IPLOT / 0 / 
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READ ( INP ) I1,I2,I3,I4,I5,I6,R1,R2,I7,I8,R3, 

6 19, 110, 111, R4, T 12, 113,114,115,116,117, 118, 119, 

6 120,121,122,123,124,125,126,127,128 

NREV » 114 
REWIND INP 
CALL USTART 
CALL UINQES (1.0, SUPRT ) 

IF ( SUPRT . NE .0.0) CALL UESCAP ( 1 . 0 , LIST, 6 . 0 , IDUM, DUM) 

C 

1000 READ ( INP, END - 2000 ) 11, 12 , 13 , 14 , 15 , 16 , Rl, R2 , 17 , 18 , R3 , 

* 19, 110, 111, R4, 112, 1 13, 114, 115, 1 16, 117, 118, 119, 

* 120,121,122,123,124,125,126,127,128 

READ ( ING, 500, END » 2000 ) IDUM1 ,GL, EL, DUM2 , DUM3 , DUM4 , 

* DUM5 , DUM6 , GUND 
500 FORMAT) 2X , 13 , 8F9 . 3 ) 

NREC * NREC + 1 

GLAT ■ 1.0E-6 * FLOAT) 17) 

ELON - 1.0E-6 * FLOAT) 18) 

HSS » 1.0E-3 * FLOAT) 113) 

IDSAT * II 
C 

IF ( 114. EQ. NREV ) THEN 

IF ( HSS. GT. EBOT. AND. HSS. LT. ETOP ) THEN 
IREC ■ IREC +• 1 
RMJD - FLOAT) 16) 

RNREV - FLOAT) I 14) 

X) IREC) - Rl 
Y(IREC) =■ HSS 
DF(IREC) » R3 
YM(IREC) » GUND 
END IF 
ELSE 

PLOT THE PREVIOUS ELEVATION PROFILE 

IF ( IPLOT . GT. MAXPLT ) GO TO 3000 
IF ( IREC. GT. 4 ) THEN . 

C 

IPLOT - IPLOT + 1 

WRITE ( OUP, 601 ) IPLOT, RNREV 

CALL PLOT ( X, Y , YM, DF , IREC , RMJD , RNREV , IDSAT , IC , NI NT , NMAX , 

6 C , F , GX , GY , DGY , DD Y , WK , XF , YF , DDT , XA , YA , DYA , IBC, 

6 SIGE , SIGT ) 

C 

END IF 

C 

NREV =* 114 
IREC * 1 

IF ( HSS. GT. EBOT. AND. HSS. LT. ETOP ) THEN 
X) IREC) » Rl 
Y(IREC) - HSS 
DF ( IREC ) = R3 
YM) IREC ) = GUND 
ELSE 

IREC = IREC - 1 
END IF 
END IF 
C 

GO TO 1000 
C 

2000 CONTINUE 
C 

IF ( IREC. GT. 4 ) THEN 

IPLOT = IPLOT + 1 

WRITE ( OUP, 601 ) IPLOT, RNREV 

CALL PLOT) X, Y, YM, DF, IREC, RMJD, RNREV, IDSAT, IC, NINT, NMAX, 

5 C , F , GX , GY , DGY , DDY , WK , XF , YF , DDT , XA , YA , DYA , IBC , 

6 SIGE, SIGT) 

END IF 

C 
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3000 WRITE ( OOP, * ) ' ' 

WRITE ( OUP, * ) ' NREC, ' 
WRITE ( OUP, * ) ' ' 

WRITE ( OUP, * ) ' I PLOT, 

C 

CALL UEND 

C 

601 TORMAT( 1 HO,//, T5, 'START PLOT 

c 

STOP 

END 

C 


RECORDS IN TILE. ' 

ELEVATION PROFILES PLOTTED. ' 


# ' ,14, ' or PASS # ' , F6 . 0 ) 


C 


SUBROUTINE PLOT 

6 

6 


x, y, 

NINT, 

DDT, 


YM, DF, IREC 
NMAX, C, F, 
XA, YA, DYA, 


, RMJD, RNREV, IDSAT, 
GX, GY, DGY, DDY, WK, 
IBC, SIGE, SIGT ) 


XC, 

XF, 


YF, 


REAL* 4 
REAL *4 
REAL *4 
REAL* 4 
INTEGER* 4 
INTEGER 
CHARACTER* 1 



IDSAT 
IBC (NMAX) 
DATE( 10) 


C ( NMAX— 1,3) 
, DDY (NMAX) 


• mtwS TH,! ABSO “ TE «"*** 


IREC 1 - IREC - 1 
REC1 m FLOAT ( IREC- 1) 

DO 110 I ■ 1, NMAX 
GX(I) - -1.0 
110 CONTINUE 


DO 111 I 
DDY ( I ) 
111 CONTINUE 
C 


- 1, IREC1 
ABS ( Y ( 1+1 ) 


Y(I) ) 


DO 112 J * 1, IREC1 

IF ( DDY ( J) . LT . SIGE ) THEN 
K ■ J 
GO TO 222 
END IF 

112 CONTINUE 
C 

222 CONTINUE 

DO 113 J ■ 1, IREC 

IF( DDY ( IREC-J) . LT. SIGE ) THEN 
IRECG = IREC-J 
GO TO 223 
END IF 

113 CONTINUE 
C 


C 


223 KB = 0 

DO 224 1=1, 

IF ( ABS ( DDY ( I 

224 CONTINUE 


IRECG- 1 
) ) -GT. SIGE) 


KB = KB 


+ 


1 


LG = 1 

GX(LG) = X ( K ) 

GY(LG) * Y ( K ) 

DGY (LG) = DF ( K ) 

888 LG = LG + 1 
K = K + 1 

IF ( ABS ( DDY ( K- 1 ) ) .GT.SIGE) THEN 
K = K + 1 

DO 225 1=1, KB- 1 

IF( ABS( DDY( K-l ) ) .GT.SIGE) K = K + 1 
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225 CONTINUE 
END IF 

GX(LG) » X(K) 

GY(LG) - Y(K) 

DGY(LG) - DF( K) 

IF ( K • LE . IRECG ) GO TO 880 
C 

C IB - 0 

C DO 234 K - 1, LG-1 

C IF(DDT(K) .GT.SIGT) THEN 

C IB - IB + 1 

C IBC(IB) * X 

C END IF 

C 234 CONTINUE 
C 

C IF(IB.EQ.l) THEN 

C KXGl * IBC(l) 

C KXG2 - LG 

C END IF 

C 

C IF(IB.EQ.l) LG - KXGl 

C 

C NA ■ 0 

C DO 999 N = 1, LG, 3 

C NA * NA + 1 

C SX * 0. 

C SY= 0 . 

C SYA - 0. 

C JC - 0 

C DO 900 L = N, N+2 

C IF(GX(L) .LT. 0. ) GO TO 246 

C JC ■ JC + 1 

C SX ■ SX + GX(L) 

C SY = SY + GY(L) 

C SYA =* SYA + DGY(L) 

C 900 CONTINUE 

C 246 RJC - FLOAT (JC) 

C XA(NA) * SX/RJC 

C YA(NA) = SY/RJC 

C DYA(NA) - SYA/RJC 

C 999 CONTINUE 
C 

C IF(IB.EQ.l) THEN 

C LG * KXG2 

C DO 998 N = KXG1+1 , LG, 3 

C NA = NA + 1 

C SX * 0. 

C SY ■ 0 . 

C SYA » 0. 

C JC = 0 

C DO 910 L * N, N+2 

C IF( GX ( L ) . LT .0.0) GO TO 357 

C JC ■ JC + 1 

C SX = SX + GX(L) 

C SY = SY + GY(L) 

C SYA = SYA + DGY(L) 

910 CONTINUE 

357 RJC * FLOAT ( JC ) 

XA(NA) = SX/RJC 
YA(NA) = SY/RJC 
DYA(NA) = SYA/RJC 
998 CONTINUE 

END IF 

....DETERMINE THE CUBIC SPLINE COEFFICIENTS C OF THE ELEVATION PROFILE 

IF ( LG . GE . 8 ) THEN 
SM = FLOAT (LG) 

CALL ICSSCU ( GX, GY , DGY , LG , SM , F , C , IC, WK, IER ) 

END IF 
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c 

c. 

c 


.START THE GRAPHIC PROCEDURES 


C 

C. 

C 


CALL USET( ' PERCEHTUNITS ' ) 

CALL USET( 'EXTRALARGE' ) 

CALL UVWPRT (0.0, 99.0, 0.0, 99.0) 
CALL UOUTLN 


ELEVATION 
: $' ) 


GEOS-3$ ’ ) 


SEASAT? ' ) 


C 

C. 

C 


C 

C 

c 


c 

c 


444 


445 


C 

C. 

C 


..PRINT TITLE, JULIAN DATE AND PASS NUMBER. 

CALL USET ( ' C JUST ' ) 

CALL USET( ' T JUST ' ) 

CALL UPSET( 'PRECISION ',6.0) 

CALL UPRINT( 50. , 97. , 

6 'ADJUSTED BLACK SEA ALTIMETER 

CALL UPRINT( 50. ,93. , 'MODIFIED JULIAN DATE 
CALL UMOVE ( 6 8 . , 9 3 . ) 

CALL UPRNT1(RMJD, 'REAL' ) 

CALL UPRINT(50.,89., 'PASS : S') 

CALL UMOVE( 58. , 89 . ) 

CALL UPRNT 1 ( RNRE V , 'REAL' ) 
ir ( IDSAT. EQ. 7502701 ) THEN 
^CALL UPRINT(50 -' 85 '' 'SATELLITE 

IF ( IDSAT. EQ. 7806401 ) THEN 
CALL UPRINT( 50 ., 85 .,’ SATELLITE 
END IF 
CALL ZTIME < DATE , 8 ) 

CALL FMOVE ( DATE ( 1 0 ) , 1 , ' $ ' ) 

CALL UPRINT( 88 . , 6 . , DATE ) 

CALL UPRINT( 88. , 3 . , 'STX/ZMAYA? ' ) 

.DRAW AND LABEL AXES (DEFAULT TIC MARKS) 

CALL USET( 'NOORIGIN' ) 

CALL UVWPRT ( 5 . , 95 . , 7 . , 85 . ) 

CALL USET( 'DSYMBOL' ) 

CALL UPSET( 'SYMBOL' ,5.0) 

CALL UPSET ( 'SZSYMBOL' , 1.0) 

CALL USET ( 'LARGE' ) 

itofvm! ^ ' , 'FRACTION OF DAY PAST MIDNIGHT?') 

J"*J*J* YLABEL ' , ' SURFACE ELEVATION IN METERS?') 

CALL USET( 'XBOTH' ) 9 ' 

CALL USET ( ' YBOTH ' ) 

CALL USET ( 'OWNSCALE' ) 

CALL UPSET ( ' TICX ' ,TCX) 

CALL UPSET ( 'TICY' ,TCY) 

XMIN - X(l) 

XMAX » X(IREC) 

YMIN * GY( 1) 

YMAX = GY( 1) 

DO 444 I = 1 , LG 

YMIN = AMIN 1 ( GY ( I ) , 

YMAX = AMAX1 ( GY( I) , 

CONTINUE 

DO 445 I = 1 , IREC 
YMIN = AMIN1 ( YM( I ) 

YMAX = AMAX1 ( YM(I) 

CONTINUE 


PROFILE? ’ ) 


YMIN 

YMAX 


YMIN 

YMAX 


555 


•DETERMINE THE INTERPOLATION INTERVALS AND PERFORM INTERPOLATION 

DX = ( GX(LG) - GX(1) ) / FLOAT) NINT) 

XF ( 1 ) = GX ( 1 ) 

DO 555 1=2, NINT+1 

XF ( I ) = XF(I-l) + DX 
CONTINUE 

IF ( LG . GE . 8 ) THEN 

CALL ICSEVU ( GX , F , LG , C , IC, XF, YF, NINT+1 , IER) 
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END IF 
C 

DO 666 I * 1 , NINT+1 

YMIN - AMIN1 ( YF ( I ) , YMIN ) 

YMAX - AMAX1 ( YF(I), YMAX ) 

666 CONTINUE 
C 

SX - 0.020 * ( XMAX - XMIN ) 

SY - 0.020 * ( YMAX - YMIN ) 

XMIN - XMIN - SX 

XMAX ■ XMAX + SX 

YMIN - YMIN - SY 

YMAX - YMAX + SY 

C SINDEX » AINT ( ABS { ALOG10 ( XMIN ) ) ) + 4. 

C XSCALE ■ 10. **S INDEX 

C XMIN - AINT (XMIN* XSCALE) /XSCALE 

C XMAX « AINT( ( XMAX * XSCALE ) +1 . ) /XSCALE 

YMIN - AINT (YMIN) - AMOD( AINT ( YMIN ) , 2 . ) 

YMAX - AINT ( YMAX ) + 2. - AMOD( AINT(YMAX), 2. ) 

CALL UWINDO< XMIN, XMAX, YMIN, YMAX) 

CALL UAXIS( XMIN, XMAX, YMIN, YMAX) 

C 

C PLOT THE DATA POINTS 

C 

DO 777 I * 1, LG 

CALL UMOVE ( GX ( I ) , GY ( I ) ) 

CALL UPRNT1 ( ' X$ ' , ' HORIZ ' ) 

777 CONTINUE 
C 

C PLOT THE SMOOTHED PROFILE 

C 

IF ( LG. GE . 9 ) THEN 
PTS - rLOAT( NINT+1 ) 

CALL USET ( ' LNULL ' ) 

CALL ULINE ( XF , YF , PTS ) 

END IF 
C 

c PLOT THE CORRESPONDING GEOID MODEL 

C 

PTS - FLOAT ( IREC ) 

CALL USET( ' DNULL ' ) 

CALL ULINE ( X, YM, PTS ) 

C 

C DISPLAY FOR SCREEN 

C 

C CALL UPAUSE 

C 

C TERMINATED THIS PLOT 

C 

CALL UERASE 
CALL URESET 
C 

RETURN 

END 

//* 

/ /*SYSLIB DD DSN=»SYS2 . IMSLD, DISP=SHR 

//* DD DSN=SYS2 . IMSLS , DISP=SHR 

/ / * 

//GO. FT03F001 DD DSN=ZMAYA . BLACKA . DATA , DISP=SHR 

//GO. FT04F001 DD DSN=ZMAYA . ALTIM . DATA ( BLKRAPP ) , DISP=SHR 

/ /* **** **** ** + * **** ★ + * ★ ***★ * it ir * * 

//* 

//* BLACK SEA DATA ADJUSTED FOR ORBIT BIAS * 

//* * 
//* *•** * * * * *■»■** **** **•* **** »»** **** »*** *•»** * * * , . 

//* 

//GO. FT06F00 1 DD SYSOUT=* 

// EXEC NOTIFYTS 
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A listing of PROGRAM SORT 


/ * J^P^^SLie 0 1 09 ' 36 ° ' 2 > ' AYAU ' TIME - <0.30), CLASS-O, MSGCLASS-X 

// EXEC rORTVC 
/SYSIN DD * 

FORMAT Or GEOS- 3 / SEASAT ALTIMETER DATA 
VARIABLE TYPE DESCRIPTION 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


14 

12 

12 

14 

14 

14 

R8 

R8 

14 

14 

R4 

14 

14 

14 

R4 

14 

14 

14 

14 

14 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 


SATELITE ID 
MEASUREMENT TYPE ( 42 
TIME SYSTEM < NM ) 

STATION NUMBER 
PREPROCESSING INDICATORS 
MODIFIED JULIAN DATE OF OBSERVATION 


OVER LAND , 43 * OVER WATER 


(GMT) 

( METERS ) 

( IE-6 DEGREES) 
( IE-6 DEGREES) 
( METERS ) 

(MM) 

(GEODYN VOL 3) 
(MM) 


FRACTION OF DAY PAST MIDNIGHT 
ALTIMETER OBSERVATION 
SATELLITE GEODETIC LATITUDE 
SATELLITE EAST LONGITUDE 
MEASUREMENT STANDARD DEVIATION 
NET INSTRUMENT CORRECTION 
METEOROLOGICAL DATA WORD 
NET MEDIA CORRECTIONS iriPl . 

GEO ID HEIGHT ABOVE REFERENCE ELLIPSOID ( METERS ) 
NET OCEAN DYNAMIC CORRECTIONS (MM) ' METERS > 
INDICATED SURFACE ELEVATION (MM) 

S/C REVOLUTION NUMBER 

ELEVATION (MARSH/MARTIN '81 (MM)) 
DOD REFERENCE RADIAL ORBIT DIFFERENCE (MM) * 

AGO 3 <CM) 

( DB ) 

WIND SPEED (CM/SEC) 

SURFACE ELEVATION PREPROCESSING WORD 
DRY TROPOSPHERIC CORRECTION (MM) 

FNOC WET TROPOSPHERIC CORRECTION (MM) 

SMMR WET TROPOSPHERIC CORRECTION (MM) 

T AM A C Dlifn T/l /v/snnMjamo a*. ’ ' 


IONOSPHERIC CORRECTION 


(MM) 

BAROTROPIC DYNAMIC SEA SURFACE CORRECTION (MM) 

ftT. TH P & D*FU T C* * ’ 


SOLID EARTH TIDE 
SCHWIDERSKI OCEAN TIDE 
PARKE OCEAN tide 


(MM) 

(MM) 

(MM) 


PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 


ETOP - 50. ) 
EBOT -10. ) 

SIGE -2. ) 

MIREC =» 150 
XMIN - 26.5 
= 42.5 
» 40.0 
= 48.0 
0.25 ) 
0.25 ) 
0.25 ) 


XMAX 
YMIN 
YMAX 
DX =* 

DY * 

CAP = 

CAP 2 = CAP* CAP ) 

NXI = ( XMAX-XMIN ) /PX + 1 ) 
NYI = ( YMAX-YMIN ) /DY + 1 ) 


INTEGER* 2 
INTEGER* 4 
REAL *8 
REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
INTEGER 


ti't4'^ 7 t£ 1 ti I1 !' I o 0 ' I21 ’ I22 ' I23 ' I24 ’ I25 ' I26 > i27 -I28 

ri’r 4,15,16,17 ' 18 ' 19 ' 110 ' 111 ' 112,113 ' 114, 115,116 

R3 R4 

GLAT, ELON, HSS, HSUM, HAVG 
GRIDX(NXI), GRIDY(NYI) 

X ( MIREC ) , Y ( MIREC ) , Z(MIREC), ER(MIREC) 

NGd”nXI?nYI? Y<MIREC) ' BZ(MIREC >' BER( MIREC), DDZ ( MIREC) 
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c 


c 


INTEGER IREC, NREC, NP, NPASS, IDSAT 
INTEGER NEG, LG, ND, NREV, MG, NTOTAL 

INTEGER INP, OUP, OUG 

DATA INP / 5 /, OUP / 6 /, OUG / 8 / 

DATA IREC / 0 /, NREC / 0 /, NEG / 0 /, ND / 0 /, NTOTAL / 0 / 

DATA NPASS / 1 / 

DATA DTR / 3 . 490658505E-2 / 

DATA HSUM / 0.0 / 


FIRST GENERATE THE GRID NET 


GRIDX(l) - XMIN 
DO 100 1-2, NXI 

GRIDX(I) - GRIDX(I-l) + DX 

100 CONTINUE 
GRIDY(l) =* YMIN 

DO 101 1-2, NYI 

GRIDY(I) - GRIDY(I-l) + DY 

101 CONTINUE 


READ ( INP 

6 

6 

IDSAT - II 
NREV - 114 
REWIND INP 


II, 12, 13, 14, 15, 16,R1,R2, 17,I8,R3, 

19, 110, II 1,R4, 112, I 13, 114, 115, 116, I 17, 118,119, 
120,121,122,123,124,125, 126,127,128 


1000 READ ( INP, END - 2000 ) I 1 , 12 , 13 , 14 , 15 , 16 , R1 , R2 , 17 , 18 , R3 , 

* 19, 110, 111, R4, I 12, I 13, 114, I 15, 116, I 17, 118, 119, 

* 120,121,122,123,124,125,126,127,128 


IF ( 114. EQ. 12793 ) GO TO 999 


NREC - NREC + 1 
GLAT - 1.0E-6 * FLOAT( 17 ) 
ELON - 1.0E-6 * FLOAT (18) 
HSS - 1.0E-3 * FLOAT ( 113 ) 


IF ( 114. EQ. NREV ) THEN 

IF ( HSS. GT. EBOT. AND. HSS. LT. ETOP ) THEN 
IREC - IREC + 1 
X(IREC) - ELON 
Y(IREC) = GLAT 
Z(IREC) - HSS 
ER( IREC) - R3 
END IF 
ELSE 


ASSIGN THE PREVIOUS ELEVATION DATA INTO AN ARRAY FOR GRIDDING 


C 

c 


IF ( IREC. GT. 4 ) THEN 

CALL SELECT ( X, Y , Z , ER, IREC , SIGE, LG , BX, BY, BZ , BER, DDZ ) 

DO 200 I = 1, LG 
ND - ND + 1 

11 = NINT( ( BX( I ) -CAP -XMIN ) /DX) 

J1 = NINT( ( BY ( I ) -CAP-YMIN ) /DY ) 

12 = NINT ( ( BX ( I ) +CAP-XMIN ) /DX ) + 2 

J2 = NINT{ ( BY ( I ) +CAP-YMIN ) /DY ) + 2 
IF ( II. LT. 1 ) II = 1 

IF ( Jl. LT. 1 ) J1 = 1 

IF ( 12 .GT. NXI ) 12 = NXI 

IF ( J2 .GT. NYI ) J2 = NYI 

DO 201 M = II, 12 
DO 202 N = Jl, J2 

R = ( GRIDX (M)-BX(I) ) * * 2 + ( GRIDY ( N ) -BY ( I ) ) * * 2 
IF ( R. LE. CAP 2 ) THEN 

NGD( M, N ) = NGD( M, N ) + 1 
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202 

201 

200 


C 

c 

c 

c 


* OUG > BX(I>,BY<I>,BZ<I), M,N, NPASS, BER(I) 

CONTINUE 

CONTINUE 

HSUM * HSUM + BZ ( I ) 

CONTINUE 

ELSE 

NPASS - NPASS - 1 
END ir 


IDSAT 
NREV = 
NPASS 
IREC ■ 


- II 
114 

■ NPASS + 1 
1 


999 


IF X(iS?; - T eloh OT ‘ AND ‘ HSS - lt - btop > THEN 
Y(IREC) - GLAT 
2 (IREC) - HSS 
ER(IREC) = R3 
ELSE . 

IREC - IREC - 1 
END IF 
END IF 

CONTINUE 
GO TO 1000 


C 

c 

c 


2000 CONTINUE 


IF ( IREC. GT. 4 ) THEN 

CALL SELECT ( X, Y, Z , ER, IREC, SIGE, LG, BX, BY, BZ , BER, DDZ ) 

DO 210 I » 1, LG 
ND - ND + 1 

11 « NINT( (BX(I)-CAP-XMIN)/DX) 

J1 » NINT ( ( BY ( I ) -CAF-YMIN ) /DY ) 

12 - NINT( (BX( I)+CAP-XMIN)/DX) + 2 
J2 - NINT< (BY(I)+CAP-YMIN)/DY) + 2 


IF ( II. LT. 
IF ( Jl. LT. 
IF ( 12 .GT. 
IF ( J2 .GT. 


1 ) 
1 ) 
NXI 
NYI 


II 
Jl - 
> 12 
) J2 


NXI 

NYI 


212 

211 

210 


DO 211 M =■ II, 12 
DO 212 N - Jl, J2 

E ~BX( I ) ) **2 + ( GRID Y (N)-BY(I) )**2 

ir ( R. LE. CAP2 ) THEN ’ " 

NGD(M,N) * NGD(M,N) + 1 

WRITE ( OUG ) BX ( I ) , BY ( I ) , BZ ( I ) , M, N, NPASS, BER(I) 


C 

C 

C 


END IF 
CONTINUE 
CONTINUE 
HSUM = HSUM 
CONTINUE 


+ BZ ( I ) 


END IF 


301 

300 


HAVG = HSUM/ FLOAT ( ND ) 

MG = -1 

DO 300 1=1, NXI 

DO 301 J = 1, NYI 

MG = MAX0 ( MG, NGD(I,J) ) 
NTOTAL = NTOTAL + NGD ( I , J ) 
CONTINUE 
CONTINUE 
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c 

c 


WRITE OUT THE DATA SET ASSOCIATED WITH EACH GRID POINT 


DO 400 I = 1, NXI 
DO 401 J - 1, NYI 
NP - NGD(I,J) 

IF ( NP. GT. 0 ) THEN 
NEG - NEG + 1 

WRITE ( OUP, 601 ) I, J, GRIDX(I), GRIDY(J), NP 
END IF 
401 CONTINUE 
400 CONTINUE 
C 

WRITE ( OUP, 604 ) 

DO 700 1*1, NXI 

WRITE ( OUP, 605 ) I, <NGD( I, J) , J=1,NYI) 

700 CONTINUE 
C 

WRITE ( OUP, 606 ) NREC, ND, NEG, MG, NTOTAL, HAVG 
C 

601 FORMAT( 1H0,T5, 'GRID < ' , 12 , ' , ' , 12 , ' ) : ' , 

4 3X, 'ELON * ' , F6 . 2 , 3X, ' GLAT - ',F6.2, 

4 5 X, 'TOTAL DATA' POINTS * ',16) 

C 602 FORMAT ( ( 1H , T5 , 2F9 . 4 , T12 . 4 ) ) 

603 FORMAT ( ( 3 ( 2F7 . 2 , F12 . 4 ) ) ) 

604 FORMAT ( 1H1,///,T20, ' DATA DISTRIBUTION MAP:',//) 

605 FORMAT ( 1H ,T5, 14 , 5X, 3313 ) 

606 FORMAT( 1H1 ,////, T5 , ’ # OF RECORDS READ: ’,18,//, 

4 T5 , ' # OF RECORDS USED IN GRIDDING: ',18,//, 

4 T5 , ' # OF NONEMPTY GRID POINTS: ',14,//, 

4 T5 , ' MAXIUM » or DATA PTS TO A GRID PT: ',15,//, 

4 T5, 'TOTAL » OF DATA RECORDS: ',110,//, 

4 T5 , ' MEAN GEOID HEIGHT OF THE BLACK SEA: ',F10.4,//, 

4 T5 , 'GRIDDED DATA HAVE BEEN WRITTEN TO BLACK4B . SORT . ' ) 

C 

STOP 
END 

SUBROUTINE SELECT ( X , Y , Z , E, IREC , SIGE, LG , GX , GY , GZ , GE , DDZ ) 
REAL *4 X(IREC), Y(IREC), Z(IREC), E(IREC) 

REAL* 4 GX(IREC), GY(IREC), GZ(IREC), GE(IREC), DDZ(IREC) 

....EVICT THE BAD DATA POINTS BY COMPARING THE ABSOLUTE DIFFERENCE 
BETWEEN ADJACENT DATA POINTS 

IRECl - IREC - 1 
REC1 = FLOAT ( IREC— 1 ) 

DO 110 1=1, IREC 

GX{ I ) = -1.0 

110 CONTINUE 
C 

DO 111 1=1, IRECl 

DDZ ( I ) = ABS ( Z(I+1) - Z(I) ) 

111 CONTINUE 
C 

DO 112 J = 1, IRECl 

IF ( DDZ ( J ) . LT . SIGE ) THEN 
K = J 
GO TO 222 
END IF 

112 CONTINUE 
C 

222 CONTINUE 

DO 113 J * 1, IREC 

IF( DDZ ( IREC-J) . LT. SIGE ) THEN 
IRECG = IREC-J 
GO TO 223 
END IF 

113 CONTINUE 
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223 KB * 0 

DO 224 1*1, IRECG-1 

IF( ABS ( DDZ ( I ) ) .GT.SIGE) 

224 CONTINUE 


KB - KB + 1 


LG - 1 
GX(LG) 
GY ( LG ) 
GZ ( LG ) 
GE(LG) 


X( K) 
Y(K) 
Z(K) 
E(K) 


225 


// 
// 
// 
II 
II 
1 1 
/ 1 
1 1 
// 

1 1* 
//* 
//* 
//* 


888 LG - LG + 1 
K - K + 1 

IF(ABS( DDZ < K-l )) .GT.SIGE) THEN 
K - K + 1 

DO 225 1*1, KB-1 

IF(ABS(DDZ(K-1) ) .GT.SIGE) K » K + 1 
CONTINUE 
END IF 

GX(LG) - X(K) 

GY (LG) * Y(K) 

GZ(LG) * Z(K) 

GE(LG) * E(K) 

IF ( K . LE . IRECG ) GO TO 888 

RETURN 

END 

EXEC LINKGOV, REGION-5000K 
*EXEC LINKGOV, REGION. LINK-5000K, REGION. GO-5000K 
* PARM* ' MAP, LIST, SIZE* ( 4096K,512K) ' 

•SYSLIB DD DSN*SYS2 . IMSLS , DISP*SHR 
DD DSN— SYS2 . IMSLD , DISP*SHR 


* ♦ * * ★ 
♦ 


* * «*♦* #4 


★ ★★* ★ ★★★ + + * + ♦* + + inn 


* 


ADJUSTED BLACK SEA DATA 


/ /S2*r r 2f r 2 01 DD OSN-ZMAYA. BLACKA. DATA, DISP-SHR 
/ /GO. FT06F00 1 DD SYSOUT** 

//GO.FT08F001 DD DSN-ZMAYA. BLACK4B . SORT, DISP-SHR 
// EXEC NOTIFYTS 
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A listing of PROGRAM WGTAVG 


//ZMAYABWA JOB ( G0109 , 360 , 2 ) , AYAU, TIME-{ 0 , 30 ) , CLASS-O, MSGCLASS-X 
/* JOBPARM LINES-60 
// EXEC FORTVC 
//SYSIN DD * 

C 

FORMAT OF GEOS-3 /SEASAT ALTIMETER DATA 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

-C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


VARIABLE 


TYPE 

14 

12 

12 

14 

14 

14 

R8 

R8 

14 

14 

R4 

14 

14 

14 

R4 

14 

14 

14 

14 

14 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 

12 


DESCRIPTION 
SATELITE ID 

MEASUREMENT TYPE ( 42* OVER LAND , 
TIME SYSTEM ( NM ) 

STATION NUMBER 

PREPROCESSING INDICATORS 

MODIFIED JULIAN DATE OF OBSERVATION 


43 


OVER WATER 


(GMT) 

( METERS ) 

( IE-6 DEGREES) 
( IE-6 DEGREES) 
( METERS ) 

(MM) 

(GEODYN VOL 3) 
(MM) 


FRACTION OF DAY PAST MIDNIGHT 
ALTIMETER OBSERVATION 
SATELLITE GEODETIC LATITUDE 
SATELLITE EAST LONGITUDE 
MEASUREMENT STANDARD DEVIATION 
NET INSTRUMENT CORRECTION 
METEOROLOGICAL DATA WORD 
NET MEDIA CORRECTIONS 

GEOID HEIGHT ABOVE REFERENCE ELLIPSOID (METERS) 

NET OCEAN DYNAMIC CORRECTIONS (MM) 

INDICATED SURFACE ELEVATION (MM) 

S/C REVOLUTION NUMBER 

MEAN SEA SURFACE ELEVATION (MARSH/MARTIN '81 (MM)) 
DOD REFERENCE RADIAL ORBIT DIFFERENCE (MM) 

H 1/3 (CM) 

AGC ( DB ) 

WIND SPEED (CM/SEC) 

SURFACE ELEVATION PREPROCESSING WORD 
DRY TROPOSPHERIC CORRECTION (MM) 

FNOC WET TROPOSPHERIC CORRECTION (MM) 

SM MR WET TROPOSPHERIC CORRECTION (MM) 


IONOSPHERIC CORRECTION 


(MM) 


BAROTROPIC DYNAMIC SEA SURFACE CORRECTION (MM) 


SOLID EARTH TIDE 
SCHWIDERSKI OCEAN TIDE 
PARKE OCEAN TIDE 


(MM) 

(MM) 

(MM) 


C 


PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 


ETOP * 50, ) 

EBOT * 10. ) 

SIGE * 2. ) 

MIREC * 150 ) 

XMIN * 26.5 ) 

XMAX =42.5 ) 

YMIN = 40.0 ) 

YMAX =48.0 ) 

DX * 0.25 ) 

DY = 0.25 ) 

MXP =100 ) 

NXI = { XMAX-XMIN ) /DX + 1 ) 

NYI = ( YMAX- YMIN ) / DY + 1 ) 


INTEGER* 2 
INTEGERM 
REAL *8 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
INTEGER 


12,13,117,118, 119,120, 121, 122, 123 f 124, 125, 126, 127,128 
11*14,15,16, 17,18,19,110, 111, 112, 113,114, 115, 116 
R 1 , R2 
R3 , R4 

GLAT, ELON, HSS, HSUM, HAVG , RSQ, W, WF, GMIN, GMAX 
EX (NXI, NYI, MXP) , GY( NXI, NYI, MXP) , H ( NXI , NYI , MXP ) 

GRIDX ( NXI ) , GRIDY(NYI) , GRID ( NXI , NYI ) 

X( MIREC), Y( MIREC), Z ( MIREC ) 

BX(MIREC), BY ( MIREC ) , BZ(MIREC), DDZ ( MIREC ) 

NGD ( NXI , NYI ) 


) 
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non 


INTEGER IREC, NREC, NP, NEG, LG, ND, NREV, MG 
INTEGER INP, OUP, OUG 
C 

DATA INP / 5 /, OUP / 6 /, OUG / 8 / 

DATA IREC / 0 /, NREC / 0 /, NEG / 0 /, ND / 0 / 

DATA HSUM / 0.0 /, GMIN / 1000.0 /, GMAX / -1000.0 / 

C 

C FIRST GENERATE THE GRID NET 

C 

GRIDX(l) - XMIN 
DO 100 1-2, NXI 

GRIDX(I) - GRIDX(I-l) + DX 

100 CONTINUE 
GRIDY(l) - YMIN 

DO 101 1-2, NYI 

GRIDY(I) » GRIDY(I-l) + DY 

101 CONTINUE 


READ ( INP 

6 

6 

NREV - 114 
REWIND INP 
C 

1000 READ ( INP, 


C 


) I1,I2,I3,I4,I5,I6,R1,R2,I7,I8,R3, 

120,121,122,123,124,125,126,127,128 


END >fl.«,I3,I4,I5, 16, R1,R2, 17,18, R3, 

19, 110, 111, R4, 112, 113, 114, 115, 116, 117, 118 119, 
120,121,122,123,124,125,126,127,128 


IF ( 114. EQ. 12793 ) 


GO TO 999 


NREC - NREC + 1 

GLAT - 1.0E-6 * FLOAT (17) 

ELON - 1.0E-6 * FLOAT (18) 

HSS - 1.0E-3 * FLOAT (113) 

IF ( 114. EQ. NREV ) THEN 

IF ( HSS. GT. EBOT. AND. HSS. LT. ETOP ) THEN 
IREC - IREC + 1 
X(IREC) - ELON 
Y( IREC) - GLAT 
Z(IREC) » HSS 
END IF 
ELSE 


....ASSIGN THE PREVIOUS ELEVATION DATA INTO AN ARRAY FOR GRIDDING 
IF ( IREC. GT. 4 ) THEN 

CALL SELECT ( X,Y,Z,IREC, SIGE, LG, BX, BY, BZ , DDZ ) 

DO 200 I - 1, LG 
ND * ND + 1 

LX - NINT ( ( BX( I ) -XMIN ) /DX ) + 1 
LY - NINT ( ( BY ( I ) -YMIN ) /DY ) + 1 
NGD( LX, LY ) = NGD( LX, LY) + 1 
EX ( LX , L Y , NGD ( LX , L Y ) ) = BX(I) 

GY(LX,LY,NGD(LX,LY) ) = BY(I) 

H ( LX , L Y , NGD ( LX , L Y ) ) = BZ ( I ) 

HSUM = HSUM + BZ ( I ) 

200 CONTINUE 

END IF 

NREV = 114 
IREC = 1 

IF ( HSS. GT . EBOT. AND. HSS. LT . ETOP ) THEN 
X ( IREC ) = ELON 
Y(IREC) = GLAT 
Z(IREC) = HSS 
ELSE 
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o o o 


1 


IREC = IREC - 
END IF 
END IF 
C 

999 CONTINUE 
GO TO 1000 
C 

2000 CONTINUE 
C 

IF < IREC. GT. 4 ) THEN 
C 

CALL SELECT ( X,Y,Z,IREC, SIGE, LG, BX, BY, BZ , DDZ ) 

C 

DO 201 I * 1, LG 
ND = ND + 1 

LX » NINT( ( BX( I ) -XMIN) /DX ) + 1 
LY » NINT( (BY( D-YMIN) /DY) + 1 
NGD( LX, LY ) - NGD( LX, LY ) + 1 
EX ( LX , LY , NGD ( LX , L Y ) ) - BX(X) 

GY ( LX , LY , NGD ( LX , L Y ) ) - BY(I) 

H ( LX , LY , NGD ( LX , LY ) ) - BZ(I) 

HSUM - HSUM + BZ<I) 

201 CONTINUE 

C 

END IF 

C 

HAVG - HSUM/FLOAT ( ND ) 

C 

MG - -1 

DO 300 I ■ 1, NXI 
DO 301 J * 1, NYI 

MG ■ MAX0 ( MG, NGD(I,J) ) 

301 CONTINUE 
300 CONTINUE 

WRITE OUT THE DATA SET ASSOCIATED WITH EACH GRID POINT 

DO 400 I - 1, NXI 
DO 401 J * 1, NYI 
NP = NGD ( I, J) 

GRID ( I , J) * -1000.0 
IF ( NP. GT. 0 ) THEN 
IF ( NP. EQ. 1 ) THEN 
GRID(I,J) = H ( I , J, 1) 

ELSE 

WF = 0.0 
W = 0.0 

DO 402 K = 1, NP 

RSQ = ( EX ( I , J , K ) -GRIDX ( I ) ) # *2 + < GY ( I , J, K ) -GRIDY ( J ) ) **2 
RSQ - 1.0/RSQ 
WF - WF + RSQ*H(I,J,K) 

W = W + RSQ 
402 CONTINUE 

GRID ( I, J) = WF/W 
END IF 

GMIN = AMIN 1 ( GMIN , GRID(I,J) ) 

GMAX = AMAX1 ( GMAX, GRID(I,J) ) 

NEG = NEG + 1 

WRITE ( OUP, 601 ) I, J, GRIDX ( I ) , GRIDY ( J ) , GRID(I,J), NP 
WRITE ( OUP, 602 ) ( EX< I , J, K ) , GY ( I , J, K ) , H ( I , J, K ) , K= 1 , NP ) 
END IF 
401 CONTINUE 
400 CONTINUE 
C 

DO 700 1=1, NXI 

WRITE ( OUG, 603 ) ( GRIDX ( I ), GRIDY ( J >, GRID ( I , J ), J= 1 , NYI ) 

700 CONTINUE 
C 

WRITE ( OUP, 604 ) 

DO 800 1=1, NXI 
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noon 


c 

c 


soo cSe 1 °“ P - 605 ’ 


WRITE ( OUP, 606 ) NREC, ND, NEG, MG, GMIN, GMAX, HAVG 

,13/) 


601 rORMAT ( 1H0 , / / , T5 , ' GRID ( ' , 12 , ' , ' , 12 , ' ) s ' 

* 3X, 'ELON - ' , 76 . 2 , 3X, 'GLAT - ',F6.2, 

- ‘ 3X,'WT. H - ' ,F10.4,5X, 'TOTAL DATA POINTS 

602 FORMAT( ( 1H , T5 , 2F9 . 4 , F12 . 4 ) ) 

603 FORMAT( ( 3 ( 2F7 . 2 , F12 . 4 ) ) ) 

604 FORMAT (1H1,///,T20,' DATA DISTRIBUTION MAP:' //) 

605 FORMAT ( 1H , T5 , 14 , 5X, 3313 ) 

606 FORMAT ( 1H1, ////,T5, 'I OF RECORDS READ: ',18,//, 

T5 , ' ♦ OF RECORDS USED IN GRIDDING: '18,//, 

T5 , ' # Or NONEMPTY GRID POINTS: ',14,//, 

T5 , ' MAXIUM I OF DATA PTS TO A GRID PT: ',13 // 

T5 , ' THE RANGE OF GEOID HEIGHT IS: ' , 2F14 . 4, // . ' 

T5 , ' MEAN GEOID HEIGHT OF THE BLACK SEA: ',Flo!4,//, 

T5, 'GRIDDED DATA HAVE BEEN WRITTEN TO BLACKA.GRID. ' ) 

STOP 
END 

SUBROUTINE SELECT ( X, Y, Z, IREC, SIGE, LG, GX, GY, GZ , DDZ ) 


REAL *4 
REAL *4 


X(IREC), Y(IREC), Z ( IREC ) 

GX(IREC), GY ( IREC) , GZ(IREC), DDZ(IREC) 


• EV ?S1,T HE BAD DATA POINTS BY COMPARING THE ABSOLUTE DIFFERENCE 
BETWEEN ADJACENT DATA POINTS 


IREC1 - IREC - 1 
REC1 - FLOAT< IREC-1) 

C 

DO 110 1*1, IREC 

GX( I ) - -1,0 

110 CONTINUE 
C 

DO 111 1*1, IREC1 

DDZ ( I ) * ABS ( Z< 1 + 1 ) 

111 CONTINUE 
C 


Z(I) ) 


DO 112 J - 1, IREC1 

IF( DDZ ( J) ,LT. SIGE) THEN 
K * J 
GO TO 222 
END IF 

112 CONTINUE 
C 

222 CONTINUE 

DO 113 J * 1, IREC 

IF ( DDZ ( IREC-J) . LT . SIGE ) THEN 
IRECG * IREC-J 
GO TO 223 
END IF 

113 CONTINUE 


223 KB * 0 

DO 224 1*1, IRECG— 1 

IF ( ABS ( DDZ ( I ) ) .GT.SIGE) 

224 CONTINUE 


KB * KB 


+ 


LG * 1 

GX(LG) = X { K ) 

GY (LG) = Y ( K ) 

GZ(LG) * Z(K) 

888 LG = LG + 1 
K = K + 1 

IT ( ABS ( DDZ { K- 1 ) ) .GT.SIGE) THEN 
K = K + I 


1 
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DO 225 1*1, KB— 1 

IF ( ABS ( DDZ ( K- 1 ) ) . GT . SIGE ) K - K + 1 
225 CONTINUE 
END IF 

GX(LG) - X{ K ) 

GY(LG) * Y(K) 

GZ ( LG ) * Z(K) 

IF { K • LE . IRECG ) GO TO 988 
C 

RETURN 

END 

//* 

// EXEC LINKGOV, REGION*5000K 
//♦SYSLIB DD DSN-SYS2 . IMSLS , DISP*SHR 

//* DD DSN*SYS2.IMSLD,DISP*SHR 

//* 

//* #**# ***★ ★ ★★★ + **% ♦ *## ♦ *** ** + ♦ ★★★★ * 
//* * 
//* ADJUSTED BLACK SEA DATA * 

//* 

//* * * ★ * **** **★* **** ***♦ * ★ ★ * * + ** * 

//* 

//GO.PT05r001 DD DSN“ZMAYA. BLACKA. DATA, DISP»SHR 
//GO. FT06F00 1 DD SYSOUT** 

/ /GO. FT08F001 DD DSN-ZMAYA. BLACKA. GRID , DISP-SHR 
// EXEC NOTIFYTS 
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non non 


A listing of PROGRAM CONTOUR 


//ZMAYABGA JOB (G0109 , 360, 10 
//TEMPLATE EXEC G38PLOT 
/ /FORT . SYS IN DD * 

C 


) , AYAU , TIME* (1,00 


CLASS-O, MSGCLASS— X 


C 


c 


c 


PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 

INTEGER 
INTEGER 
INTEGER 
INTEGER 
REAL *4 
REALM 
REALM 
REAL* 4 
CHARACTER* 1 


BIAS - 0.0 ) 

XMIN - 26.5 ) 

XMAX - 42.5 ) 

YMIN - 40.0 ) 

YMAX - 48.0 ) 

DX - 0.25 ) 

DY - 0.25 ) 

CHGH - 50. ) 

CLOW » 10. ) 

CINC - 1. ) 

CLAB - 2. ) 

NLON - ( ( XMAX-XMIN ) /DX + 1 ) ) 

NLAT - ( ( YMAX- YMIN ) /DY + 1 ) ) 

INP, NRU, OUP 
NPROJ, LINEH, LINEV 
LIST ( 6 ) 

MASK (NLON, NLAT) 

LON (NLON), LAT(NLAT), GRID ( NLON, NLAT) 
WKCNTR ( NLON , NLAT ) , XP ( NLON ) , YP ( NLAT ) 
PENS ( 7 ) 

PLAT , PLONG, VLAT1 , VLAT2 , VLONG1 , VLONG2 
DATE (10) 


DATA 

DATA 

DATA 

DATA 

DATA 

DATA 


INP / 5 /, OUP / 6 /, NRU / 0 / 
GMAX / -999.0 /, GMIN / 1000.0 / 
PENS / 7*4.0 / 

LIST / 8, 110, 75, 0, 7, 10 / 
NPROJ / 16 / 

PLAT / 90. /, PLONG / 0. / 


DO 100 I - 1, NLON 
READ ( INP, 510, END 
100 CONTINUE 
C 


777 ) ( LON ( I ) , LAT ( J) , GRID ( I , J ) , J— 1 , NLAT ) 


777 CONTINUE 
C 


DO 111 I - 1, NLON 
DO 112 J - 1, NLAT 

IF ( GRID( I , J ) . NE. -1000.0 ) THEN 
NRU - NRU + 1 

GRID ( I , J ) - GRID(I,J) -BIAS 
GMAX - AMAX1 ( GRID(I,J), GMAX ) 
GMIN - AMIN 1 ( GRID ( I , J ) , GMIN ) 
END IF ' 

112 CONTINUE 
111 CONTINUE 


WRITE ( OUP, 610 ) GMIN, GMAX, NRU 

WRITE ( OUP ,620) ( ( LON ( IX ) , LAT( IY ) , GRID( IX, IY ) , IY=1 , NLAT ) , IX— l , NLON ) 


START GRAPHIC PROCEDURES 


CALL USTART 


CALL UINQES (1.0, SUPRT ) 

IF ( SUPRT. NE . 0 . 0 ) CALL UESCAP(1.0 
CALL USET( ' PERCENTUNITS ' ) 

CALL UVWPRT( 0. 0, 99 . 0, 0 . 0, 99 . 0 ) 

CALL UOUTLN 


, LIST, 6 . 0 


IDUM , DUM ) 


PUT ON LABELS 
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CALL USET ( ' C JUST ' ) 

CALL USET( ' TJUST ' ) 

CALL USET ( 'YELLOW' ) 

CALL USET ( 'SOFTWARE' ) 

CALL USET( 'EXTRALARGE' ) 

CALL UPSET ( 'VERTICAL SIZE' ,2.0) 

CALL UPSET( 'HORIZONT SIZE ',1.6) 

CALL UFONT( 'TROM' ) 

CALL UPRINT(50. ,96. , 

fc 'WT. AVG. GEOIDAL HEIGHT OF THE BLACK SEA$ ' ) 

VLATl * YMIN 

VLAT2 ■ YMAX 

VLONGl * XMIN 

VLONG2 * XMAX 

LINEH * IFIX( YMAX- YMIN) 

LINEV - IFIX(XMAX-XMIN) 

CALL WDMAP ( NPRO J , PLAT , PLONG , VLAT 1 , VLAT2 , VLONG 1 , VLONG2 , 
fc LINEH, LINEV) 

LABEL AXIS AND WRITE LON AND LATS ON PLOT 

DLONG - ABS ( VLONG2 - VLONGl) / LINEV 

LINES - LINEV + 1 

DIX * 100.0 / LINEV 

YLONGI - -3.5 

DO 420 1*1, LINES, 2 

V * VLONGl + (DLONG * (I - 1)) 

XLATIT * DIX * (I - 1) 

CALL UMOVE( XLATIT, YLONGI) 

CALL UPRNT1 ( V, ' REAL ' ) 

420 CONTINUE 

WRITE LABEL 

CALL USET( ' BJUS ' ) 

CALL UMOVE ( 5 0 . , -7 . 5 ) 

CALL UPRNT1 ( ' LONGITUDES ' , ' HORIZ ' ) 

DLAT * ABS ( VLAT2 - VLATl) / LINEH 
LINES * LINEH + 1 
DIY * 100.0 / LINEH 
XLT * -3.5 

DO 430 1*1, LINES, 2 

V = VLATl + (DLAT * (I - 1)) 

YLO * DIY * ( I - 1 ) 

CALL UMOVE (XLT, YLO) 

CALL UPRNT1(V, 'REAL' ) 

430 CONTINUE 

WRITE LABEL 

CALL USET ( ' MJUS ' ) 

CALL UMOVE (-7. 5, 5 0.0) 

CALL UPRNT1( 'LATITUDES' , 'VERTI' ) 

....WRITE DATE AND ID ON PLOT 

CALL ZTIME ( DATE , 8 ) 

CALL FMOVE ( DATE ( 1 0 ) , 1 , ' $ ' ) 

CALL UPRINT ( 7 8 . , -8 . , DATE ) 

CALL UPRINT (95.,— 8.,' STX/ ZMAYAS ' ) 

....START CONTOURING PROCEDURES 

CALL USET ( 'NOMI ' ) 

CALL USET ( 'SMOOTH' ) 

CALL UPSET ( 'CLOWEST' , CLOW) 

CALL UPSET ( 'CINCREMENT' ,CINC) 

CALL UPSET ( 'CHIGHEST' , CHGH ) 
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CALL UPSET ( 'CLABEL' , CLAB) 

C 

CALL USET( ' NOAXES ’ ) 

CALL USET( ' NOXLABEL ' ) 

CALL USET( 'NOYLABEL' ) 

c CONVERT LAT/LON TO PLOT COORDS. 

C 

DO 222 I-l,NLON 

CALL WDPOS (1.0, LON ( I ) , XP ( I ) , DUM ) 

222 CONTINUE 

DO 223 J » 1, NLAT 

« 3 coSSoS DPOS 1 “ T 1 J 1 • 1 • 0 ■ 00H ' ” 1 J " 
c 

C WRITE(OUP, 630) 

C WRITE (OUP, 620) ( <XP(IX),YP<IY),GRID(IX,IY>,IY-l,NLAT),IX-l,NLON) 


C 

c 

c 

c 


FX - FLOAT (NLON) 

FY » FLOAT (NLAT) 

CALL UVWPRTfO. 0,99. 0,0. 0,99.0) 

CALL USET ( 'OWNSCALE' ) 

CALL UPSET ( 'EXCEPTION' ,-1000.0) 

CALL USET ( ' DNULL ' ) 

CALL UPCNTR ( GRID , XP , YP , WKCNTR, FX , FY , PENS ) 

CALL UPAUSE 
CALL UEND 


F12.4, 


510 FORMAT ( ( 3 ( 2F7 . 2 , F12 . 4 ) ) ) 

61 ° £ r ? R ^ T | 1 S!h T 5' ' THE RANGB OF GRID VALUES IS FROM 
* IN UNMASKED RECORDS.'// 

4 T10 , LONGITUDE, LATITUDE AND SURTACE HEIGHT' /) 

620 FORMAT ( ( 1H , 4 ( 2F7 . 2 , 2X, F12 . 4 , IX) )) ,/) 

c 630 rORMAT <lHl,T15, 'COORDINATES AND GRID VALUES FOR CONTOURING',/) 

STOP 

END 

//* 

//GO. FT05F001 DD DSN-ZMAYA. BLACKA. GRID, DISP*SHR 

//* **** **** .**. *.** .*** .*** 

//* * **** **** * 

//* ADJUSTED BLACK SEA GRID DATA USING WT. AVERAGING * 

//. * .... * 

//GO. FT06F001 DD SYSOUT=* 

//GO.TT19F001 DD DSN=SYS2 . WRLDATA2 , DISP*SHR, LABEL=»( , , , IN) 

// EXEC NOTIFYTS 
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A listing of PROGRAM GEOID 


//ZMAYABSH JOB ( GO 109 , 36 0 , 2 ) , AYAU, TIME-< 2 , 30 ) , CLASS-O, MSGCLASS-X 
// EXEC FORTVC 
//SYSIN DD * 

C 

PARAMETER ( MPASS >142 ) 

PARAMETER ( NREFP * 73 ) 

PARAMETER ( KUTOFF * 3 ) 

PARAMETER ( LCTV * 1 ) 

PARAMETER ( MTR » 160 ) 

PARAMETER ( XMIN - 26.5 ) 

PARAMETER ( XMAX = 42.5 ) 

PARAMETER ( YMXN » 40.0 ) 

PARAMETER < YMAX - 48.0 ) 

PARAMETER ( DX * 0.25 ) 

PARAMETER ( DY - 0.25 ) 

PARAMETER ( LCT * 33 ) 

PARAMETER ( NXI - ( XMAX-XMIN ) /DX + 1 ) 

PARAMETER ( NYI = ( YMAX-YMIN ) /DY + 1 ) 

PARAMETER ( LEN - ( MTR* ( MTR+1 ) ) /2 ) 

PARAMETER ( MSYM * ( MPASS* ( MPASS- 1 )) /2 ) 

C 


C 


REAL *4 
REAL *4 
REAL *4 
REAL *8 
REAL *8 
REAL *8 
REAL* 8 
REAL *8 
REAL'S 
REAL* 8 
INTEGER 
INTEGER 


SE ( MTR ) 

GRIDX(NXI), GRIDY(NYI) 

GRID (NXI, NYI ) , VARM( NXI , NYI ) 

RGRIDX(NXI), RGRIOY ( NYI ) 

SY(MTR), CY(MTR) 

EX ( MTR ) , GY ( MTR ) , GH(MTR), SGY(MTR), CGY(MTR) 
CM(LEN), CMX(LEN), CMS ( LEN ) , DM ( LEN ) , DLL ( LEN ) 
CP ( MTR ) , SRC ( MTR) , AP(MTR) 

SD(LCT), CTV ( LCT , 3 ) , CV(LCT), SC(LCT-l f 3) 

DTR, ELON, GLAT, XI, Yl, SY1, CY1, PVAL, VAR 
ID I AG ( MTR ) , NPASS < MTR ) 

INP, OUP, OUG, OUH, OUC, OUE 


DATA INP / 5 /, OUP / 6 /, OUG / 8 /, OUH / 9 / 
DATA OUC / 10 /, OUE / 11 / 

DATA NIER / 0 / 

DATA DTR / 3 . 490658504D-2 / 

DATA GMIN / 1000.0 /, GMAX / -999.0 / 

DATA VMIN / 1000.0 /, VMAX / -999.0 / 


GENERATE A GRID NET 

GRIDX(l) » XMIN 
DO 100 I > 2, NXI 

GRIDX(I) * GRIDX(I-l) + DX 
RGRIDX(I) > DBLE ( GRIDX ( I ) ) * DTR 

100 CONTINUE 
GRIDY(l) = YMIN 

DO 101 I =• 2, NYI 

GRIDY(I) = GRIDY(I-l) + DY 
RGRIDY ( I ) = DBLE ( GRIDY ( I ) ) * DTR 

SY ( I ) * DSIN ( RGRIDY ( I ) ) 

CY ( I ) = DCOS ( RGRIDY ( I ) ) 

101 CONTINUE 

DO 102 1*1, NXI 

DO 103 J = 1, NYI 
GRID ( I , J ) = -1000.0 
VARM( I , J ) = -1000.0 
103 CONTINUE 

102 CONTINUE 

....READ IN THE COVARIANCE TABLE 
DO 110 1=1, LCT 

READ ( INP, 501 ) SD ( I ) , ( CTV ( I , J ) , J= 1 , 3 ) 
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CV(I) . 
110 CONTINUE 


CTV ( I , LCTV ) 


£ DETERMINE A CUBIC SPLINE COEFFICIENT MATRIX SC 

ISC - LCT - 1 

CALL ICSCCU ( SD, CV, LCT, SC, ISC, IER ) 

£ WRITE OUT THE VALUES OF THE COVARIANCE TABLE USED 

WRITE ( OUP, 601 ) 

c WRITE ( OUP, 602 ) (SD(I),CTV(I, LCTV), 1-1, LCT) 

WRITE ( OUP, 600 ) 

C 

£ READ IN THE ERROR COVARIANCE MATRIX FROM PASS BIAS ADJUSTMENT 

REAO < OUE , 502 ) ( DM ( M ) , M* 1 , MPASS 1 ) 

£ CONSTRUCT THE COVARIANCE MATRIX DIAGONAL ADDRESS 

DO 120 1*1, MTR 

IDIAG(I) - (I*(I+l))/2 
120 CONTINUE 
C 

c READ IN THE SORTED DATA 

C 

READ < OUG , Rl, R2 , R3, II, 12, 13, R4 
IKEEP ■ II 
JKEEP - 12 
REWIND OUG 

0,000 SSS i : IT > **« «> «• «. i3. R4 

GLAT » DBLE ( R2 ) * DTR 
HSS - R3 


C 

c 


c 

c 


IF ( II. EQ. IKEEP. AND. 12. EQ. JKEEP ) 

IREC * IREC + 1 
EX(IREC) - ELON 
GY (IREC) - GLAT 
GH(IREC) ■ DBLE (HSS) 

NPASS ( IREC ) * 13 
SE(IREC) - R4 

IF ( NPASS (IREC). EQ. NREFP ) 

IF ( NPASS ( IREC ) . GT. NREFP ) 


THEN 


NPASS ( IREC) 
NPASS ( IREC) 


NPASS ( IREC) 


ELSE 


678 


IF ( IREC. GT. KUTOFF ) THEN 
XI * RGRIDX( IKEEP) 

Y1 - RGRIDY( JKEEP) 

SY1 » SY( JKEEP) 

CY1 * CY( JKEEP) 

JREC = ( IREC*( IREC+1 ) ) /2 
JC - JC + 1 

IF ( JC.GE.571 .AND. JC.LE.590 ) THEN 
WRITE ( OUP, * ) 

WRITE ( OUP, 678 ) ( EX( LK ), GY ( LK ), GH ( LK ), NPASS ( LK ), SE ( LK ) , 

LK= 1 , IREC ) 

FORMAT) ( 1H ,3F12. 4, I5,F12.4>> 

WRITE ( OUP, * ) 

END IF 

CALL PREDIC ( PVAL.VAR, EX, GY, GH, SGY, CGY, IREC 

NPASS, SE, DM, DLL, XI , Y1 , SY1 , CY1 , 

CM, CMS, CMI , JREC, IDIAG, CP, SRC , AP, 

SD , CTV , CV , SC , LCT , ISC , 

„„ LCTV, DTR, IER, NIER ) 

IF ( IER. EQ. 0 ) THEN 
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GRID) I KEEP, JKEEP) ■ PVAL 
VARM ( I KEEP , JKEEP ) - VAR 
ELSE 

JC - JC - 1 

WRITE ( OUP, * ) ' IER * ' , IER, ' PVAL =',PVAL, 

* ' VAR - ' , VAR 

END IF 

WRITE ( OUP, 603 ) IKEEP, JKEEP, GRIDX( IKEEP) ,GRIDY( JKEEP) , 

* GRID ( IKEEP, JKEEP) , VARM ( IKEEP, JKEEP) , IREC,JC 
END IF 

C 

IKEEP - II 
JKEEP - 12 
IREC * 1 
EX(IREC) = ELON 
GY ( IREC ) - GLAT 
GH(IREC) ■ DBLE(HSS) 

NPASS ( IREC ) - 13 
SE(IREC) - R4 

IF ( NPASS (IREC). EQ. NREFP ) NPASS (IREC) » 0 
IF ( NPASS (IREC). GT. NREFP ) NPASS (IREC) =» NPASS (IREC) - 1 
C 

END IF 
C 

GO TO 1000 

c 

2000 CONTINUE 

c 


IF ( IREC. GT. KUTOFF ) THEN 
XI * RGRIDX(Il) 

Y 1 » RGRIDY ( 12 ) 

SY1 * SY( JKEEP) 

CY1 - CY( JKEEP) 

JREC = ( IREC*( IREC+1) ) /2 
JC - JC + 1 

CALL PREDIC ( PVAL, VAR, EX, GY , GH , SGY, CGY, IREC, 

* NPASS, SE, DM, DLL, XI , Y1 , SY1 , CY1 , 

* CM, CMS, CMI, JREC, IDIAG, CP, SRC, AP, 

* SD,CTV,CV,SC,LCT, ISC, 

* LCTV, DTR, IER, NIER ) 

IF ( IER. EQ. 0 ) THEN 

GRID( 11,12) * PVAL 

VARM( 11,12) - VAR 

ELSE 

JC ■ JC - 1 

WRITE ( OUP, * ) 'IER = ' , IER, ' PVAL =',PVAL, 

* ' VAR = ' , VAR 

END IF 

WRITE < OUP, 603 ) 11,12, GRIDX( II), GRIDY ( 12 ) , 

* GRID( I 1 , 12 ) , VARM( I 1 , 12 ) , IREC,JC 
END IF 


....DETERMINE THE RANGE OF THE PREDICTED VALUES AND THEIR VARIANCE 

DO 300 1=1, NXI 

DO 301 J = 1, NYI 

IF ( GRID ( I , J ) . NE. -1000.0 ) THEN 
GMIN = AMIN 1 ( GMIN , GRID) I, J) ) 

GMAX = AMAX1 ( GMAX, GRID) I, J) ) 

END IF 

IF ( VARM( I , J) . NE. -1000.0 ) THEN 
VMIN = AMIN 1 ( VMIN , VARM) I, J) ) 

VMAX = AMAXl ( VMAX , VARM) I, J) ) 

END IF 
301 CONTINUE 
300 CONTINUE 


PRINT OUT THE GRIDDED VALUES AND THE VARIANCE 


NXY 


NXI *NYI 
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WRITE ( OUP, 604 ) 

DO 400 I - 1, NXI 
WRITE ( OUP, 605 ) 
WRITE ( OUH, 606 ) 
400 CONTINUE 


NIER, JC, NXY, GHIN, GMAX 

( GRIDX( I ) , GRIDY ( J) , GRID( I , J) , J*1 , NYI ) 
( GRIDX ( I ) , GRIDY ( J ) ,GRID( 1 , J) \ J-iInYI) 


WRITE ( OUP, 607 ) VMIN, VMAX 
DO 410 I - 1, NXI 
WRITE ( OUP, 605 ) 

WRITE ( OUC, 606 ) 

410 CONTINUE 


( GRIDX ( I ) , GRIDY ( J) , VARM( I, J) , J=1 , NYI ) 

<GRIDX(I),GRIDY(J,,VARM(I^J),J-1,NYI» 


501 FORMAT ( 

502 FORMAT ( 
500 FORMAT ( 

601 FORMAT ( 

602 FORMAT ( 

603 FORMAT ( 


604 FORMAT ( 


605 FORMAT ( 

606 FORMAT ( 

607 rORMAT( 


IX, D10 . 4 , 3D20 . 13 ) 

(4D20. 13)) 

1H1 ) 

1H1,T5, 'VALUES OF THE COVARIANCE TABLE USED:’ //) 

<1H , 5X,D10.4,5X,D20.13) ) 

( Z \' ' ’ ‘ ,X2 ' ’ > ' ' 3x » 'ELON - ’ , F6 . 2 , 2X, 
/ r «-2,5X, 'GEOID HEIGHT » ',F10.4,' +- ' F8 2 5X 
RECORD •: ', 15 , 5X, ' COUNTER: ’,15,/) ' 2,5X ' 

1H1,T5, 14, ' GRID POINTS HAVING NEGATIVE VARIANCE' // 

T5, 'THERE ARE ',14,' POINTS GRIDDED OUT OF A NET' 

' OF ’,15,' GRID POINTS. ’ ,//,T8, 'THE RANGE OF THE ' 

( 1H ,5<2F?T^.^f? HT IN METERS 
<3(2F7.2,F12.4) ) ) 

1H1,T5,'THE RANGE OF THE VARIANCE IN CM IS : ' , 2F12 . 4 , // ) 


STOP 

END 

SUBROUTINE PREDIC ( PVAL , VAR , EX , GY , GH, SGY, CGY, IREC, 

NPA SS, SE, DM, DLL, XI , Y1 , SYl , CY 1 , 
CM, CMS , CMI , JREC , IDIAG, CP , SRC , AP , 
SD , CTV , CV , SC , LCT , ISC , 

LCTV , DTR , IER, NIER ) 


REAL *4 
REAL *8 
REAL *8 
REAL* 8 
REAL *8 
REAL *8 
REAL* 8 
INTEGER 


SE(IREC) 

EX(IREC), GY(IREC), GH(IREC), SGY (IREC), CGY f IREC! 
CP(IREC), SRC(IREC), AP(IREci ' ’ ' CGY < IREC > 

CMI ( JRE C), CMS (JREC), DM(JREC), DLL (JREC) 
SD(LCT), CTV ( LCT , 3 ) , CV(LCT), SC(LCT-1,3) 

Dl, D2 , CSD , DSD, VDSD, AVG, DTR 
XI, Yl, SYl, CY 1 , PVAL, VAR 
IDIAG ( IREC) , NPASS ( IREC) 


C 

C. 

C 


.RESET THE ERROR COUNTER 
IER - 0 


C 

C. 

c 


.CLEAR THE COVARIANCE MATRIX ARRAY 

CALL DCLR( CM, JREC ) 

CALL DCLR( CMI, JREC ) 

CALL DCLR( CMS, JREC ) 

CALL DCLR( DLL, JREC ) 

.DETERMINE THE AVERAGE VALUE OF THE DATA SET 


AVG =0.0 

DO 100 1=1, IREC 

AVG = AVG + GH ( I ) 

100 CONTINUE 

avg = avg / d Float ( irec ) 


* DIAG0NAL TERMS OF THE SIGNAL COVARIANCE MATRIX 
REMOVE THE AVERAGE VALUE FROM DATA SET TO MAKE IT ZERO MF AN 
AND SET UP THE DATA VECTOR IN DOUBLE PRECISION 

DO 110 1=1, IREC 

CM ( IDIAG ( I ) ) = CTV { 1 , LCTV ) + DBLE ( SE ( I ) *SE ( I ) ) 
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GH(I) = GH(I) - AVG 
110 CONTINUE 
C 

C EVALUATE ALL THE SINE AND COSINE TERMS 

C 

DO 200 I - 1, IREC 
SGY(I) * DSIN ( GY ( I ) ) 

CGY(I) = DCOS ( GY ( I ) ) 

200 CONTINUE 
C 

c CONSTRUCT THE OFT- DIAGONAL TERMS OF THE SIGNAL COVARIANCE MATRIX 

C 

DO 300 J - 2, IREC 
DO 301 I » 1, J-l 
C 

CSD - SGY ( I ) *SGY ( J) + CGY ( I ) *CGY< J) *DCOS < EX< I ) -EX( J ) ) 

DSD - DACOS(CSD) /DTR 

CALL ICSEVU ( SD, CV, LCT, SC, ISC, DSD, VDSD, 1, IER ) 

K - ( J*< J-l) )/2 + I 
CM( K ) = VDSD 
C 


301 CONTINUE 
300 CONTINUE 
C 

C ASSEMBLE THE ERROR COVARIANCE MATRIX FROM BIAS ADJUSTMENT 


C 

DO 310 J-l, IREC 
DO 311 I » 1, J 
C 

K * ( J*( J-l) )/2 + I 

IF ( NPASS(I). NE. 0. AND. NPASS(J). NE. 0 ) THEN 
Ml - MIN0 ( NPASS(I), NPASS(J) ) 

M2 ” MAX0 ( NPASS(I), NPASS(J) ) 

L - (M2* (M2-1 ) )/2 + Ml 
DLL ( K ) = DM(L) 

ELSE 

DLL ( K ) * 0.0DO 
END IF 
C 


311 CONTINUE 
310 CONTINUE 
C 

C ASSEMBLE THE COMPLETE COVARIANCE MATRIX AND 

C DUPLICATE THE COVARIANCE MATRIX FOR STORAGE 


C 

DO 320 I - 1, JREC 

CM( I ) - CM( I ) + DLL ( I ) 
CMS ( I ) - CM( I ) 

320 CONTINUE 


C 

C INVERT THE COVARIANCE MATRIX 

C 

CALL LINV1P ( CMS, IREC, CMI, IDGT, D1,D2, IER ) 

C CALL SOLVE ( CMS, SRC, IDIAG, IREC, .TRUE., .FALSE. ) 

C CALL SYMINV ( CMS, SRC, IDIAG, IREC, CMI ) 

C 

IF ( IER. EQ. 0 ) THEN 
C 

c CONSTRUCT THE COVARIANCE VECTOR 

C 

DO 330 1=1, IREC 

C 


CSD = SY 1 *SGY ( I ) + CY1*CGY( I ) *DCOS( X1-EX( I ) ) 

DSD = DACOS(CSD) /DTR 

CALL ICSEVU ( SD, CV, LCT, SC, ISC, DSD, VDSD, 1, IER ) 
CP(I) = VDSD 
C 


330 CONTINUE 
C 

C DETERMINE THE PREDICTED VALUE AT A GRID POINT 
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c 

c 

c 

c. 

c 


c 

c 

c 

c 


c 

c 


c 

c 

c 


c 

c 


CALL VMULSF { CMI, IREC, GH,1,IREC, SRC,IREC ) 

CALL DCLR ( SRC, IREC ) 

CALL VSMXVT ( CMI , GH, SRC, IDIAG, IREC ) 

PVAL - DOT ( CP, SRC, IREC ) + AVG 

DETERMINE THE VARIANCE OF THE PREDICTED VALUE 

CALL VMULSF ( CMI, IREC, CP, 1, IREC, AP, IREC ) 

CALL VMULSF ( CM, IREC, AP,1,IREC, SRC, IREC ) 

CALL DCLR ( AP, IREC ) 

CALL DCLR ( SRC, IREC ) 

CALL VSMXVT ( CMI, CP, AP, IDIAG, IREC ) 

CALL VSMXVT ( CM, AP, SRC, IDIAG, IREC ) 

VAR - CTV{ 1 , LCTV ) - 2 . ODO*DOT( AP, CP , IREC ) + DOT ( AP , SRC , IREC ) 

IF ( VAR. GT. 0.0 ) THEN 
VAR * 100. 0D0*DSQRT( VAR) 

ELSE 

IER - 1 

NIER » NIER + 1 
END IF 

END IF 

RETURN 

END 

FUNCTION DOT ( A, B , N ) 

Compute the dot product of the two N-vectors A and B. 


INTEGER 
REAL* 4 
REAL* 8 

DOT*0 . 0 


N, I 
DOT 

A( 1) ,B< 1) 


C 

C 


DO 100 I-1,N 

100 DOT*DOT + A ( I ) *B ( I ) 

RETURN 

END 

SUBROUTINE SYMINV ( A, B, IDIAG, NEQ, AINV ) 
IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A ( 1 ) , B(l), AINV(l), IDIAG< 1 ) 

K * 0 

DO 100 1*1, NEQ 

DO 101 L * 1, NEQ 
B(L) = 0.0D0 

101 CONTINUE 

B ( I ) - 1.0D0 

® OLVE _< A ' B ' IDIAG, NEQ, .FALSE., .TRUE. ) 
DO X 1 0 J — 1 , X 
K * K + 1 
AINV(K) = B(J) 

110 CONTINUE 
100 CONTINUE 


RETURN 

END 

SUBROUTINE SOLVE ( A, B , IDIAG, NEQ, FACT, BACK) 

* U factorization of the symmetric matrix 
stored m A, if FACT = TRUE; and solve A * X = B if BACK = TRUE. 


co ™P« c ted-column form of the upper triangular 
part of the coefficient matrix. After factorization, it 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c. 


c 


c 

c 

c 


c 

c 

c. 

c 


c 

c 

c 

c 

c 

c 


contains D and U. 

B Right-hand-aide vector . After backsubstitution, it 

contains the solution. 

IOIAG Addresses of the diagonal terms in A. 

NEQ Number of equations 

FACT If FACT « TRUE, then factor A; otherwise do not factor A. 

BACK If BACK ■ TRUE, reduce B and backsubstitute; otherwise 

do not solve the equations. 


IMPLICIT REAL* 8 (A-H,0-Z) 

LOGICAL FACT, BACK 
DIMENSION A( I) ,B< 1) , IDIAG( 1) 


Factor A, reduce B 
JR » 0 

DO 400 J =» 1, NEQ 
JD - IDIAG<J) 

JH - JD - JR 
IS - J - JH + 2 

IF (JH .LT. 2) GOTO 390 

IF (FACT) THEN 

IF (JH .GT. 2) THEN 

. Reduce column J rows IS to J— 1: do not divide by row diagonal 

K « JR + 2 
ID - IDIAG( IS - 1) 

DO 100 I * IS, J-l 
IR - ID 
ID - IDIAG(I) 

IH - MIN ( ID-IR-1, I-IS+1) 

IF (IH.GT.0) A(K) * A(K) - DOT( A( K-IH ) , A( ID-IH ) , IH ) 

K =» K + 1 

100 CONTINUE 

ENDIF 

. Divide by row diagonal, and reduce diagonal term in column J 

IR * JR + 1 
K * J - JD 

DO 200 I » IR, JD- 1 
ID * IDIAG ( K+I ) 

IF (A(ID).EQ.O.O) GOTO 200 
D =» -A ( I ) 

A( I ) * A( I ) /A( ID ) 

A( JD) * A( JD) + D*A( I ) 

200 CONTINUE 

ENDIF 

. Reduce RHS 

IF (BACK) B ( J ) * B ( J ) - DOT ( A < JR+ 1) ,B(IS-1) ,JH-1) 

390 JR * JD 
400 CONTINUE 


IF (.NOT. BACK) RETURN 


Divide by diagonal pivots 
DO 700 I = 1 , NEQ 
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ID - IDIAG(I) 

IF <A(ID) .NE.O.O) 
700 CONTINUE 
C 

C Backsubstitute 
C 


B(I) 


B(I)/A(ID) 


J - NEQ 

JD ■ IDIAG(J) 

801 D - -B(J) 

J « J - 1 

ir (J.LE.O) RETURN 


JR - IDIAG(J) 

IF (JD-JR.GT.l) THEN 
IS * J - JD + JR + 2 
K - JR - IS + 1 
DO 810 I » is, J 
810 B( I ) =» B(I) + A( I+K ) *D 

ENDIF 
C 

JD - JR 
GOTO 801 
C 

c 

END 

SUBROUTINE VSMXVT ( A, B, 

C THIS SUBROUTINE FORMS C » 

C STORED IN PROFILE FORM, 

C DIAGONALS IN A. 

C 


C, JDIAG, NEQ ) 

C + A*B WHERE A 
B, C ARE VECTORS 


x nr Lav.ii u \ n,v— 

DIMENSION A<1), 0(1), C<1), JDIAG(l) 


100 


200 


JS - 1 
DO 200 
JD - 
IF ( 
0J 


J - 1, NEQ 
JDIAG ( J) 

JS. LE. JD ) THEN 
B< J) 

AB * A( JD) *BJ 
IF ( JS. NE. JD ) THEN 
JB » J - JD 
JE - JD - 1 
DO 100 JJ ■ JS, JE 

AB * AB + A( JJ) *B< JJ+JB) 

C ( JJ+JB ) * C ( JJ+JB ) + A( JJ ) *BJ 
CONTINUE 
END IF 

C(J) - C(J) + AB 
END IF 
JS ■ JD + 1 
CONTINUE 


RETURN 

END 

SUBROUTINE DCLR ( A , N A ) 

INTEGER NA, I 
REAL *8 A ( 1 ) 

C 

DO 100 1*1, NA 

A< I) = 0. 0D0 
100 CONTINUE 
C 


RETURN 

END 

//* 

// EXEC LINKGOV, REGION=*5 0 0 OK 

// + 
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IS A SYMMETRIC MATRIX 
, AND JDIAG LOCATES THE 



//SYSLIB DD DSN-SYS2 . IMSLD, DISP-SHR 
//* DD DSN-SYS2 . IMSLS, DISP-SHR 

//* 

//GO . FT05F00 1 DD DSN-ZMAYA. ALTIM. DATA< BEl 80 ), DISP-SHR 

//* 

//* **** **** * *** **** **** **** **** * 

//* 

//* BLACK SEA DATA AFTER SORT * 

//* * 
/ /* ★ *** #★** **** **** **** **** * 

//* 

/ /GO . FT06F00 1 DD SYSOUT-* 

/ /GO . FT08F00 1 DD DSN-ZMAYA. BLACK4 . SORT, DISP-SHR 
/ /GO. FT09F001 DD DSN-ZMAYA. BLACKC4H . GRID, DISP-SHR 
//GO.FTlOrOOl DD DSN-ZMAYA. BLACKC4H . COVM, DISP-SHR 
//GO.FT11F001 DD DSN-ZMAYA. BLACK. PCOV, DISP-SHR 
// EXEC NOTIFYTS 
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A listing of PROGRAM EMPCOV 


t 


/ /ZMAYAB36 JOB 
// EXEC FORTVC 
/ /SYSIN DD * 

C 

C 


(G0109 , 360,2), AYAU , TIME- (00,30), CLASS-F , MSGCLASS-X 


C 


PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 
PARAMETER ( 


XMIN - 26.5 ) 

XMAX - 42.5 ) 

YMIN - 40.0 ) 

YMAX - 48.0 ) 

DXY - 0.25 ) 

M - ( XMAX- XMIN ) /DXY + 1 
N - (YMAX- YMIN) /DXY + 1 
MN - N ) 


) 

) 


C 


IMPLICIT REAL* 8 (A-Z) 

DIMENSION G(M,N), DG(M,N), H(M, N) 
D( 0 : MN) 

CHH ( 0 : MN ) , CHHNS ( 0 : MN ) , 
CHG ( 0 : MN ) , CHGNS ( 0 : MN ) , 
CGG ( 0 : MN ) , CGGNS ( 0 : MN ) , 
INP, OUP, IMH, INH, IMG 


DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

INTEGER 


, DH(M,N) 

CHHEW ( 0 : MN ) 
CHGEW(0:MN) 
CGGEW ( 0 : MN ) 
, ING, OUF 


DATA INP / 5 /, 
* IMG / 9 /, 

C 
C 

c READ IN THE MODEL 

C 

DO 100 J - N, 1, 
READ ( IMH, * ) 
READ ( IMH, 501 
100 CONTINUE 
C 


OUP / 6 /, 
ING / 10 / 


IMH / 7 /, INH / 8 /, 
OUF / 11 / 


GEO I DAL HEIGHT AND GRAVITY ANOMALY 


-1 

RL 


) ( H ( I , J ) , 1—1 , M) 


C 

C 

c 


DO 101 J » N, 1, 
READ ( IMG, * ) 
READ ( IMG, 501 
101 CONTINUE 


-1 

RL 

) (G(I,j),i=*i,m) 


READ IN THE GEOIDAL HEIGHT AND GRAVITY ANOMALY 


DATA 


C 


C 

C 

C 


C 

C 

c 


DO 110 I - 1, M 
READ ( INH, 502 
110 CONTINUE 


(DH(I,J;,J-1,N) 


DO 111 I - 1, M 
READ ( ING, 502 
111 CONTINUE 


(DG( I, J) , J-1,N) 


REMOVE THE MODEL VALUE FROM THE DATA 


DO 120 I = 1, M 
DO 121 J = 1, 
IF ( DH ( I , J ) 
IF ( DG( I , J) 
121 CONTINUE 
120 CONTINUE 


N 

• NE. 
.NE. 


- 1000.0 ) 
- 1000.0 ) 


DH( I, J) 
DG( I, J) 


PERFORM THE CONVOLUTION 


DO 200 K =* 0, MN - 1 


SHH = 0. 

SHG - 0. 

SGG = 0. 

CHHNS ( K ) = 0. 
CHGNS ( K ) = 0. 
CGGNS ( K ) = 0 . 


DH( I, J) - H( I, J) 
DG< I, J) - G< I, J) 
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DO 210 J « 1, N 
DO 211 I - K+l, M 

IF ( DH(I-K, J) .NE. -1000.0 .AND. DH( I , J) . NE. -1000 . 0 ) 
SHH - SHH + DH< I-K, J) *DH< I, J) 

END IF 

IF ( DH( I-K, J) .NE. -1000. 0 .AND. DG( I , J) . NE. - 1000 . 0 ) 
SHG - SHG + DH< I-K, J) *DG( I, J) 

END ir 

IF ( DG(I-K,J) .NE. -1000.0 .AND. DG( I, J) .NE. -1000 . 0 ) 
SGG - SGG + DG( I-K, J) *DG( I , J) 

END IF 
211 CONTINUE 

210 CONTINUE 

Cl * ( (M-K) *N) 

CHHNS(K) = SHH / Cl 
CHGNS(K) = SHG / Cl 
CGGNS(K) » SGG / Cl 

SHH - 0. 

SHG - 0. 

SGG - 0. 

CHHEW(K) = 0. 

CHGEW(K) ■ 0. 

CGGEW(K) - 0. 

DO 220 I - 1, M 
DO 221 J - K+l, N 

IF ( DH< I, J-K) .NE.-1000. 0 .AND. DH ( I , J ) . NE . -1000 . 0 ) 
SHH » SHH + DH< I, J-K) *DH< I, J) 

END IF 

IF ( DH< I, J-K) .NE. -1000.0 .AND. DG( I , J) . NE . - 1000 . 0 ) 
SHG * SHG + DH ( I , J-K ) +DG( I , J) 

END ir 

IF ( DG(I, J-K). NE. -1000.0 .AND. DG< I, J) . NE. -1000 . 0 ) 
SGG = SGG + DG( I, J-K) +DG( I, J) 

END IF 
221 CONTINUE 

220 CONTINUE 

C2 - (M*(N-K) ) 

CHHEW(K) - SHH / C2 

CHGEW(K) = SHG / C2 

CGGEW(K) - SGG / C2 


D(K) * K * DXY 
CHH(K) * ( Cl+CHHNS ( K ) 
CHG(K) - ( C 1+CHGNS ( K ) 
CGG(K) = ( Cl*CGGNS(K) 
C 

200 CONTINUE 
C 

WRITE ( OUP, 601 ) 

WRITE ( OUP, 602 ) ( D( 

WRITE ( OUF , 603 ) ( D( 


+ 

C2*CHHEW(K) 

) 

/ 

( 

Cl 

+ 

C2 

) 

+ 

C2*CHGEW< K) 

) 

/ 

( 

Cl 

+ 

C2 

) 

+ 

C2*CGGEW( K) 

) 

/ 

( 

Cl 

+ 

C2 

) 


), CHH(I), CHG(I), CGG(I), I - 0, 
), CHH(I), CHG(I), CGG(I), 1=0, 


C 


501 FORMAT ( 

502 FORMAT ( 
601 FORMAT ( 


602 FORMAT ( 

603 FORMAT ( 


( 6 ( 2X, F10 . 3 ) ) ) 

( 3 ( 14X, F12 . 4 ) ) ) 

1H1,///,T5, 'SPHERICAL DISTANCE ', T30 C ( NN ) IN M**2' 
T48 , ' C( NG ) IN M*MGAL ' ,T67, 'C(GG) IN MGAL * * 2 


T5 , ' ' , T30 , ' 

T48, ' ' , T6 7 , ' — 


( 1H ,T15,F5.2,T30,F12. 4,T50,F12.4,T7 0,F12.4) ) 
( IX, D10 . 4, 3D2 0 . 13 ) ) 


STOP 

END 

//* 

// EXEC LINKGOV, REGION=5000K 

//* 

//GO. FT06F00 1 DD SYSOUT=* 

/ / GO .FT07F001 DD DSN=ZMAYA. ALTIM. DATA( BRH36 ) , DISP=SHR 


THEN 

THEN 

THEN 


THEN 

THEN 

THEN 


MN-1 ) 
MN-1 ) 


/, 

,/) 
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//GO.FT08F001 DD DSN— ZMAYA. ALTIM. DATA( RDH300 ) 
//GO.FT09r001 DD DSN-ZMAYA. ALTIM. DATA/ BRG36 ) , 
//GO. FT10F001 DD DSN-ZMAYA. ALTIM. DATA< RDG300 ) 
//GO.FTIlFOOl DD DSN-ZMAYA. ALTIM. DATA( BEI036 ) 

// EXEC NOTITYTS 


, DISP-SHR 
DISP-SHR 
, DISP-SHR 
, DISP-SHR 
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A listing of PROGRAM GRAVEN 


/ /ZMAYABOR JOB ( GO 109 ,360,2) , AYAU , TIME“( 04,00) , CLASS»0, MSGCLASS«X 
/♦JOBPARM LINES-60 
// EXEC FORTVC 
//SYSIN DO * 

C 

PARAMETER ( MODEL ■ 180 ) 

PARAMETER ( ANOISE - 0.25 ) 

PARAMETER ( ITM - 0 ) 

PARAMETER ( CAP - 1.00 ) 

PARAMETER ( IPRT - 0 ) 

PARAMETER ( CNVRGE * 1 . 0E-4 ) 

PARAMETER ( SCX = 1.00 ) 

PARAMETER ( KUTOFF ■ 3 ) 

PARAMETER ( MTR - 600 ) 

PARAMETER ( XMIN -26.5 ) 

PARAMETER ( XMAX - 42.5 ) 

PARAMETER ( YMIN » 40.0 ) 

PARAMETER ( YMAX - 48.0 ) 

PARAMETER ( XGl - 26.5 ) 

PARAMETER ( XG2 -42.5 ) 

PARAMETER ( YG1 - 40.0 ) 

PARAMETER ( YG2 - 48.0 ) 

PARAMETER ( DX • 0.25 ) 

PARAMETER ( DY » 0.25 ) 

PARAMETER ( NXI - <XG2-XGl)/DX + 1 ) 

PARAMETER ( NY I - ( YG2-YG1 ) /DY + 1 ) 

PARAMETER ( Ml - ( XMIN-XG1 ) /DX + 1 ) 

PARAMETER ( M2 - ( XMAX-XG1 ) /DX + 1 ) 

PARAMETER ( N1 - ( YMIN-YG1 ) /DY + 1 ) 

PARAMETER ( N2 ■ ( YMAX-YG1 ) /DY + 1 ) 

PARAMETER ( LCT * NY I ) 

PARAMETER ( INC - 4. * CAP + 1 ) 

PARAMETER ( CAP2 » CAP*CAP ) 

PARAMETER ( LEN - ( MTR* ( MTR+1 ) ) /2 ) 

PARAMETER ( LCTF • (2*LCT-1) ) 

PARAMETER ( NFFT - ( 6 *LCTF/2 ) + 150 ) 


REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
REAL* 4 
REAL* 4 
REAL*4 
REAL* 4 
REAL *4 
REAL *8 
REAL *8 
REAL *8 
REAL *8 
REAL *8 
REAL* 8 
REAL *8 
REAL *8 
REAL *8 
REAL *8 
REAL *8 
REAL* 8 
REAL *8 
REAL *8 
REAL *8 
REAL *8 
COMPLEX* 16 


GRIDX(NXI), GRIDY(NYI) 

GUI ( NXI , NYI ) , HVARI ( NXI , NYI ) 

GU ( NXI , NYI ) , HVAR( NXI , NYI ) 

GA< NXI , NYI ) , GVAR( NXI , NYI ) 

DIFF ( NXI , NYI ) 

GRA300 ( NXI , NYI ) , GUN300 ( NXI , NYI ) 

GRA180 ( NXI , NYI ) , GUN 1 80 ( NXI , NYI ) 

GRA036 ( NXI , NYI ) , GUN036 ( NXI , NYI ) 

GRAN (NXI, NYI ) , GUND< NXI , NYI ) 

WKCG ( NXI , NYI ) , WKCH ( NXI , NYI ) 

CHHNS ( 0 : NYI- 1 ) , CHHEW( 0 : NYI- 1 ) 

CHGNS ( 0 : NYI- 1 ) , CHGEW( 0 : NYI-1 ) 

CGGNS ( 0 : NYI- 1 ) , CGGEW( 0 : NYI-1 ) 

RGRIDX ( NXI ) , RGRIDY ( NYI ) 

RMSH ( 0 : ITM) , RMSG(0:ITM) 

RMSH1( 0: ITM) , RMSG1(0:ITM) 

RMSH2< 0: ITM) , RMSG2(0:ITM) 

RMSH3( 0 : ITM) , RMSG3(0:ITM) 

SY(MTR), CY(MTR), SE(MTR) 

EX (MTR), GY (MTR), GH(MTR), SGY(MTR), CGY ( MTR ) 
CM(LEN), CMI(LEN), CMS ( LEN ) 

CP (MTR), SRC (MTR), AP(MTR) 

SD(LCT), CTV ( LCT , 3 ) 

CVUU(LCT), SCUU ( LCT- 1,3) 

CVUG(LCT), SCUG ( LCT- 1,3) 

CVGG(LCT), SCGG ( LCT- 1 , 3 ) 

WKF(NFFT), A(LCTF) 

PCVUU(LCTF), PCVUG(LCTF) , PCVGG(LCTF) 

DTR, ELON, GLAT, XI, Yl, SY1, CY1, PVAL, VAR 
X(LCTF) 
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o o o n 


C 


c 

c 

c 


c 

c 

c 


INTEGER 

INTEGER 

INTEGER 


IDIAG(MTR), MASK ( NXI , NYI ) 
INP, OUP, INH, INV, INU1 , 
OU G, OVG, OUH, OVH* OUD, 


, IWKF(NFFT) 
INGl , INU2 , 
OUB 


ING2 , 


INM 


DATA 

DATA 

DATA 

DATA 

DATA 

DATA 

DATA 

DATA 

DATA 

DATA 

DATA 

DATA 


HE/, 5 .O' 3 *'*/'™* / 8 /, INV / 9 / 


INGl 

ING3 


11 /, 
15 / 


/ 18 /, OUH / 19 /, OVH / 20 /, 


INU1 / 10 /, 

INU3 / 14 /, 

INM / 16 / 

OUG / 17 /, OVG 
OUB / 22 / 

NIER / 0 / 

DTR / 3 . 4906585040-2 / 

GMING / 1000.0 /, GMAXG / -999.0 

VMING / 1000.0 /, VMAXG / -999.0 

GMINH / 1000.0 /, GMAXH / -999.0 

VMINH / 1000.0 /, VMAXH / -999.0 


INU2 / 12 /, ING2 / 13 / 


OUD 


READ IN THE BLACK SEA MASK 


READ ( INM, 505 ) 
WRITE ( OUP, 612 ) 
WRITE ( OUP, 505 ) 


( <MASK(I,J),I-M1,M2),J-N2,N1,-1) 

( ( MASK( I, J) ,I»M1,M2) , J»N2,N1,-1) 


GENERATE A GRID NET Or SINE AND COSINE 


100 


101 


'■* < A ; — AVj» A 

RGRIDX( I ) « DBLE ( GRIDX(l) 
DO 100 I » 2, NXI 

GRIDX(I) - GRIDX ( I — 1 ) + D 
RGRIDX(I) - DBLE ( GRIDX ( 
CONTINUE 
GRIDY(l) » YGl 
RGRIDY(l) m DBLE ( GRIDY<1) 
SY(1) * DSIN ( RGRIDY ( 1 ) ) 

CY ( 1 ) - DCOS ( RGRIDY ( 1 ) ) 

DO 101 I - 2, NYI 

GRIDY(I) - GRIDY(I-l) + 
RGRIDY ( I ) - DBLE ( GRIDY ( ] 
SY ( I ) - DSIN ( RGRIDY ( I ) ) 

CY ( I ) - DCOS ( RGRIDY < I ) ) 
CONTINUE 


* DTR 


) * DTR 
* DTR 

) * DTR 


’ *** ^HE MODEL GEOIDAL HEIGHT AND GRAVITY 

THE OSU 180X180, 36X36 AND 300X300 MODELS 


ANOMALY 


102 

C 


DO 102 J « NYI , 1, 
READ ( INU 1 , * ) 
READ ( INU 1 , 503 ) 
CONTINUE 


-1 

RL 

( GUN 180( I # J) , 1*1, NXI 


) 


103 

C 


DO 103 J « NYI, 1, 
READ ( INGl, * ) 
READ ( INGl, 503 ) 
CONTINUE 


-1 

RL 

( GRAl 8 0 ( I, J) , 1*1, NXI 


) 


104 

C 


DO 104 J » NYI, 1, 
READ ( INU2, * ) 
READ ( INU2, 503 ) 
CONTINUE 


-1 

RL 

( GUN036 ( I, J) , 1=1 , NXI 


) 


105 

C 


DO 105 J » NYI, 1, 
READ ( ING2, * ) 
READ ( ING2 , 503 
CONTINUE 


-1 

RL 

(GRA036 ( 1 , J) , 1=1 , NXI ) 


C 


106 


DO 106 J * NYI, 1, 
READ { INU3, * ) 
READ ( INU3, 503 ) 
CONTINUE 


-1 

RL 

{ GUN3 00(1, J) ,1 = 1, NXI 


) 


21 / 
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DO 107 J a NYI, 1, -1 
READ ( ING3 , * ) RL 

READ ( ING3, 503 ) (GRA300( I, J) , 1*1, NXI) 

107 CONTINUE 
C 

C READ IN THE GRIDDED GEOID UNDULATION DATA AND VARIANCES 

C 

DO 108 I -Ml, M2 

READ ( INH, 504 ) (GRIDX( I) , GRIDY( J) , GUI( I , J) , J*N1 , N2 ) 

108 CONTINUE 
C 

DO 109 I * Ml, M2 

READ ( INV, 504 ) ( GRIDX( I ) , GRIDY ( J) , HVARI ( I , J) , J*N1 , N2 ) 

109 CONTINUE 
C 

IF ( IPRT .NE. 0 ) THEN 
WRITE ( OUP, 608 ) 

WRITE ( OUP, 609 ) ( ( GRIDX ( I ) , GRIDY ( J) , GRA036 ( I , J) , GUN036 ( I , J ) , 

* GRA1 80 ( I , J) , GUN 1 80 ( I , J) , GRA3 00 ( I , J ) , 

* GUN300 ( I , J) , GUI ( I , J ) , HVARI ( I , J ) , 

* J*1 , NYI ) ,1*1, NXI ) 

END IF 

C 

c SELECT THE APPROPRIATE OSU MODEL 

C 

IF ( MODEL . EQ. 180 ) THEN 
C 

DO 110 I * Ml, M2 
DO 111 J * Nl, N2 

GRAN ( I , J ) * GRA180 ( I , J ) 

GUND ( I , J ) * GUN 1 8 0 ( I , J ) 

111 CONTINUE 

110 CONTINUE 
C 

ELSE IF ( MODEL .EQ. 36 ) THEN 
C 

DO 112 I - Ml, M2 
DO 113 J - Nl, N2 

GRAN ( I , J ) * GRA036 ( I , J) 

GUND ( I , J ) * GUN036 ( I , J ) 

113 CONTINUE 

112 CONTINUE 
C 

ELSE IF ( MODEL .EQ. 300 ) THEN 
C 

DO 114 I = Ml, M2 
DO 115 J * Nl, N2 

GRAN ( I , J ) » GRA3 0 0 ( I , J ) 

GUND( I , J ) * GUN 3 00 ( I , J ) 

115 CONTINUE 

114 CONTINUE 
C 

ELSE 

C 

WRITE ( OUP, 999 ) MODEL 
STOP 
C 

END IF 
C 

C READ IN THE COVARIANCE TABLE 

C 

READ ( INP, 501 ) ( SD ( I ) , ( CTV ( I , J ) , J* 1,3), 1=1, LCT) 

c 

c SCALE THE CORELATION LENGTH OF THE COVARIANCE FUNCTION 

C 

DO 149 1=1, LCT 

SD<I) = SD ( I ) * DBLE(SCX) 

149 CONTINUE 
C 

c CONSTRUCT THE COVARIANCE MATRIX DIAGONAL ADDRESS 
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c 

c 

c. 

c 


DO 116 I i 
IDIAG ( I ) 
116 CONTINUE 


1, MTR 

• <IMI + l))/2 


C 

C. 

C 


•START AN ITERATION PROCESS TO REFINE THE GRAVITY ANOMALY 
DO 900 IJK - 0, ITM 


150 


C 

C. 

c 


c 

c. 

c 


c 

c. 

c 

c 


DO 150 I 
CVUU(I) 
CVUG< I) 
CVGG( I) 
CONTINUE 


1, LCT 
CTV (1,1) 
CTV(I,2) 
CTV (1,3) 


•DETERMINE THE POWER SPECTRA OF THE COVARIANCE FUNCTIONS 

CALL POWSPC ( A, CVUU, X, PCVUU, LCT, LCTF, NFFT IWKF WKF \ 

CALL POWSPC ( A, CVUG, X, PCVUG LCT LCTF NFFT IWKf' WKF 

CALL POWSPC ( A, CVGG, X, PCVGG; LC^ IZtY, NfS,' WKF ) 

.DETERMINE A CUBIC SPLINE COEFFICIENT MATRIX SC 

* CSCCU < SD, CVUU, LCT, SCUU, LCT-1, IER ) 

* CSCCU < SD » CVUG, LCT, SCUG, LCT-1, IER ) 

CALL ICSCCU ( SD, CVGG, LCT, SCGG, LCT-1, IER ) 

.WRITE OUT THE VALUES OF THE COVARIANCE FUNCTIONS AND THE 
CORRESPONDING POWER SPECTRA AND THE 


WRITE ( OUB, 631 ) 

WRITE ( OUB; 501 ) 

WRITE ( OUB, 632 ) 

WRITE ( OUB, 501 ) 

WRITE ( OUP, 601 ) 

WRITE ( OUP, 602 ) 

WRITE ( OUP, 634 ) 

WRITE ( OUP, 602 ) 


IJK 

(SD(I),CTV<I,1),CTV<I,2),CTV(I,3),I-1,LCT) 

IJK 

( SD( I ) , PCVUU ( I ) , PCVUG ( I ) , PCVGG < I ) , 1 = 1 , LCT) 
UK 

IJK ( I ) , CTV( 1,1), CTV( 1,2), CTV (1,3), 1=1 , LCT ) 
( SD( I ) , PCVUU ( I ) , PCVUG ( I ) , PCVGG ( I ) , 1=1 , LCT ) 
.TRANSFORM THE GEOIDAL HEIGHT TO GRAVITY ANOMALY 


WRITE ( OUP, 600 ) 

WRITE ( OUP, ♦ ) 

WRITE ( OUP, * ) 


GEOIDAL HEIGHT TO GRAVITY ANOMALY ; IJK 


RMS ^ 
RMS 1 
RMS 2 
RMS3 
JCG 


0. 

- 0. 
* 0. 
■ 0, 
0 


DO 200 I * Ml, M2 
DO 201 J - Nl, N2 

GA ( I , J ) = -1000.0 


GVAR 

< I 

,J) 

= -1 

000. 

IF ( 

MASK ( 

I, J) 

.NE 

GUI { 

I, J) 

.NE 

• -n 

11 

3 

I - 

INC 


12 

3 

I + 

INC 


Jl 

3 

J - 

INC 


J2 

= 

J + 

INC 


IF 

( 

11. 

LT. 

1 ) 

IF 

( 

12 . 

GT. 

NXI 

IF 

( 

Jl. 

LT. 

1 ) 

IF 

( 

J2 . 

GT. 

NYI 


.AND. GRAN ( I , J ) . NE . -1000.0 .AND. 


THEN 


II = 1 
) 12 = NXI 

Ji = 1 
) J2 = NYI 
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IREC 


0 



c 


c 


c 

211 

210 

C 


c 


c 

201 

200 

C 


c 

c 

c 

c 


c 


DO 210 M - II, 12 
DO 211 N ■ Jl, J2 


IF (GUI(M,N) .NE. -1000. .AND. GUND(M,N) .NE. -1000. ) THEN 
RS - (GRIDX(I)-GRIDX(M) )**2 + (GRIDY( J ) -GRIDY ( N ) ) **2 
IF ( RS. LT. CAP2 ) THEN 
IREC - IREC + 1 
EX(IREC) - RGRIDX(M) 

GY (IREC) » RGRIDY(N) 

GH(IREC) ■ DBLE ( GUI(M,N) - GUND ( M, N ) ) 

SE(IREC) - DBLE ( 0 . 0 1*HVARI ( M, N ) + ANOISE ) 

END IF 
END IF 


CONTINUE 

CONTINUE 


IF ( IREC. GT. KUTOFF ) THEN 
XI » RGRIDX(I) 

Y 1 * RGRIDY(J) 

SY1 =■ SY(J) 

CY1 * CY(J) 

JREC - ( IREC*( IREC+1) ) /2 
JCG - JCG + 1 

CALL PREDIC ( 1, PVAL,VAR, EX,GY,GH, SGY, CGY, IREC, 

SE, XI, Y1,SY1,CY1, 

CM,CMS,CMI, JREC, IDIAG, CP, SRC , AP , 




LCT, DTR, IER, NIER ) 
IF ( IER. EQ. 0 ) THEN 

GA( I , J ) - PVAL + GRAN ( I , J) 

GVAR( I , J ) - VAR 
RMS * RMS + PVAL* PVAL 
D1 - GA( I , J ) - GRA1 90 ( I , J) 

GA( I, J) - GRA036 ( I , J) 

G A ( I , J ) - GRA300(I,J) 

“ RMS1 + D1*D1 
D2*D2 
D3*D3 


RMS 2 
RMS 3 


D2 - 
D3 =■ 

RMS 1 
RMS 2 
RMS 3 

ELSE 

JCG = JCG - 1 
IF ( IPRT .NE. 

' PVAL 

END IF 

IF ( IPRT .NE. 0 


0 ) WRITE ( OUP, * 


■' , PVAL, ' 


VAR 


) 'IER =' , IER, 
' , VAR 


END IF 


) WRITE ( OUP, 603 ) I , J , GRIDX ( I ) , 
GRIDY ( J ) , GA ( I , J ) , GVAR( I , J) , PVAL, IREC, JCG 


END IF 


CONTINUE 

CONTINUE 

RMSG(IJK) = SQRT ( RMS / FLOAT (JCG) ) 

RMSGl(IJK) * SQRT ( RMS1 / FLOAT (JCG) ) 

RMSG2 ( I JK ) = SQRT ( RMS 2 / FLOAT (JCG) ) 

RMSG3 ( IJK ) = SQRT ( RMS 3 / FLOAT (JCG) ) 

WRITE ( OUP, 614 ) RMSG( IJK) ,RMSG3< IJK) ,RMSG1( IJK) , RMSG2 < IJK) 

TRANSFORM BACK FROM GRAVITY ANOMALY TO GEOID HEIGHT 
WRITE ( OUP, 600 ) 

WRITE ( OUP, * ) ' GRAVITY ANOMALY BACK TO GEOIDAL HEIGHT: ' 

X JK ' 

WRITE ( OUP , * ) 

RMS = 0. 

RMS 1 = 0. 
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c 


RMS 2 - 0. 

RMS 3 * 0. 

JCH - 0 

DO 250 I ■ Ml, M2 
DO 251 J * N 1 t N2 


C 

c 


c 

c 

c 


c 

261 

260 

c 


GU(I,J) - -1000.0 
HVAR(I,J) - -1000.0 


IF ( MASK(I,J) .NE. 0 .AND. GUND(I,J) 
G A ( I , J ) .NE. -1000.0 .AND. GRAN(I,J) 

11 - I - INC 

12 - I + INC 

Jl ■ J - INC 


J2 

» 

J ♦ 

INC 





15 

( 

11. 

LT. 

i ) 

11 - 

1 


IF 

< 

12. 

GT. 

NXI 

) 12 

m 

NXI 

15 

( 

Jl. 

LT. 

1 > 

Jl « 

1 


IF 

( 

J2. 

GT. 

NYI 

) J2 

s 

NYI 


NE. -1000.0 .AND. 
.NE. -1000.0 ) THEN 


IREC - 0 


DO 260 M - II, 12 
DO 261 N - Jl, J2 


.NE.-1000. . AND. GRAN (M,N) . NE.-1000. ) T 
RS ■ (GRIDX( I ) -GRIDX( M) ) *»2 + ( GRIDY ( J ) -GRIDY ( N 
ir { RS. LT. CAP2 ) THEN 
IREC - 


EX (IREC) 
GY ( IREC ) 
GH( IREC) 
SE( IREC) 
END IF 
END IF 


LT. CAP2 
IREC + 1 


THEN 
) >**2 


RGRIDX(M) 

- RGRIDY(N) 

- DBLE ( GA ( M , N ) - GRAN ( M , N ) 

- DBLE ( GVAR(M,N) ) 


CONTINUE 

CONTINUE 


Xi - RGRIDX(I) 

Y 1 * RGRIDY ( J) 

SY1 » SY ( J) 

CY1 » CY ( J ) 

JREC - ( IREC* ( IREC+1 ) ) /2 

JCH - JCH + 1 

CALL PREDIC ( 2, PVAL,VAR 


EQ. 


IF ( IER. 

GU( I, J) 

HVAR( I, J) 

RMS - 

D 1 = 

D2 =■ 

D3 = 

RMS 1 
RMS2 
RMS 3 
ELSE 

JCH = JCH - 1 
IF ( IPRT .NE. 

' PVAL = ' 

END IF 

IF ( IPRT .NE. 0 


rvftL^AK, EX , GY , GH , SGY , CG Y , I REC 
SE, XI, Y1,SY1,CY1, 

CM, CMS , CMI , JREC , IDIAG , CP , SRC , AP , 
SD , CTV , CVUU , SCUU , CVUG , SCUG , CVGG . S 
LCT, DTR, IER, NIER ) 


LCT, DTR, IER, NIER ) 
0 ) THEN 
PVAL + GUND ( I , J ) 

100. *VAR 
RMS + PVAL* PVAL 
GU ( I , J ) - G"N 1 80 ( I , J ) 

GU( I, J) - GUN036 ( I, J) 

GU( I , J ) - GUN 3 0 0 ( I , J ) 

D 1 *D1 
D2*D2 
D3*D3 


SCGG, 


RMS 1 
RMS 2 
RMS 3 


0 ) WRITE 
, PVAL , ' 


OUP, 
VAR = 


* ) 

; ' , VAR 


IER =' , IER, 


END IF 


GRIDY ( J ) , 


WRITE ( 
GU ( I , J ) , 


OUP, 617 ) 
HVAR ( I , J ) , 


I , J , GRIDX ( I ) , 
PVAL, IREC, JCH 
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noon n n noon non non 


c 


END IF 


C 

251 CONTINUE 

250 CONTINUE 

RMSH(IJK) - SQRT ( RMS / FLOAT (JCH) ) 

RMSHl(IJK) - SQRT ( RMS1 / FLOAT ( JCH) ) 

RMSH2 ( IJK) - SQRT ( RMS 2 / FLOAT (JCH) ) 

RMSH3( IJK) - SQRT ( RMS 3 / FLOAT (JCH) ) 

WRITE ( OUP, 614 ) RMSH ( IJK ) , RMSH3 ( IJK ) , RMSH1 ( IJK ) , RMSH2 ( IJK ) 

CHECK CONVERGENCE CRITERION 

IF ( RMSG(IJK) .LT. CNVRGE ) GO TO 777 

REDETERMINE THE EMPIRICAL COVARIANCE FUNCTION 

CALL EMPCVF ( SD,CTV, NXI , NYI , DX, NYI , GRAN , GA , GUND,GUI, 

* WKCG, WKCH , CHHNS , CHHEW, CHGNS , CHGEW , CGGNS , CGGEW ) 

USE THE RESULTANT GEOIDAL HEIGHT AND GRAVITY ANOMALY AS A 

WORKING MODEL FOR THE NEXT ITERATION STEP 

DO 270 1*1, NXI 

DO 271 J - 1, NYI 
GUND( I, J) » GU( I, J) 

271 CONTINUE 

270 CONTINUE 

DO 280 1*1, NXI 

DO 281 J * 1, NYI 
GRAN ( I , J ) * GA ( I , J ) 

281 CONTINUE 

280 CONTINUE 

900 CONTINUE 


....WRITE OUT THE SUMMARY OF THE ITERATION. 

777 IF ( IJK .GT. ITM ) IJK * ITM 

WRITE ( OUP, 618 ) MODEL 

WRITE ( OUP, 621 ) ( L, RMSG(L), RMSH(L), RMSG3 ( L ) , RMS H 3 ( L ) 

RMSGl(L), RMSHl(L), RMSG2 ( L ) , RMSH2 ( L ) , 

L * 0, IJK ) 

....DETERMINE THE RANGE OF THE PREDICTED VALUES AND THEIR VARIANCE 

DO 300 I - Ml, M2 
DO 301 J = Nl, N2 

IF ( GA ( I , J ) . NE. -1000.0 ) THEN 
GMING = AMIN 1 ( GMING, GA(I,J) ) 

GMAXG = AMAXl ( GMAXG, GA(I,J) ) 

END IF 

IF ( GVAR( I , J ) . NE. -1000.0 ) THEN 
VMING = AMIN 1 ( VMING , GVAR(I,J) ) 

VMAXG = AMAXl ( VMAXG, GVAR(I,J) ) 

END IF 

IF ( GU ( I , J ) . NE. -1000.0 ) THEN 
GMINH = AMIN 1 ( GMINH , GU ( I , J ) ) 

GMAXH = AMAXl ( GMAXH , GU(I,J) ) 

END IF 

IF ( HVAR ( I , J ) . NE. -1000.0 ) THEN 
VMINH = AMIN 1 ( VMINH , HVAR(I,J) ) 

VMAXH = AMAXl ( VMAXH, HVAR(I,J) ) 

END IF 
301 CONTINUE 
300 CONTINUE 
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noon 


C PRINT 

C 

NXY - 
WRITE ( 
DO 400 
WRITE 
WRITE 

400 CONTINUE 


OUT THE GRIDDED VALUES AND THE VARIANCE 
( M2 -M 1+1 ) * ( N2-N1+1 ) 


OUP, 604 ) 
I - Ml, M2 
( OUP, 605 
( OUG, 606 


NIER, JCG, NXY, GMING, GMAXG 

) (GRIDX( I) ,GRIDY( J) ,GA( I, J) , J-N1,N2 ) 
) ( GRIDX( I ) , GRIDY ( J) , GA( I , J) , J=N1 , N2 ) 


WRITE ( OUP, 615 ) 
DO 401 I - Ml, M2 
WRITE < OUP, 605 
C WRITE ( OUH, 606 

401 CONTINUE 
C 


NIER, JCH, NXY, GMINH, GMAXH 

) ( GRIDX ( I ) , GRIDY ( J ) , GU ( I , J ) , J«N1,N2) 

) ( GRIDX ( I ) , GRIDY ( J ) , GU ( I , J ) , J»N 1 , N2 ) 


WRITE ( OUP, 607 ) 
DO 410 I - Ml, M2 
WRITE ( OUP, 605 
WRITE ( OVG, 606 
410 CONTINUE 


VMING, VMAXG 

) { GRIDX ( I) , GRIDY (J) , G VAR ( I , J ) , J=N1,N2 ) 

) ( GRIDX ( I ) , GRIDY ( J) , GVAR( I , J) , J»Nl , N2 ) 


WRITE ( OUP, 616 ) 
DO 411 I ■ Ml, M2 
WRITE ( OUP, 605 
C WRITE ( OVH, 606 

411 CONTINUE 


VMINH, VMAXH 


) ( GRIDX ( I) , GRIDY 

) ( GRIDX ( I ) , GRIDY 


( J) ,HVAR(I, J) , J-N1,N2 
( J) , HVAR< I, J) , J«N1,N2 


) 

) 


DETERMINE THE DIFFERENCE BETWEEN 
AND THE MOOEL VALUES USED 


THE PREDICTED GRAVITY ANOMALIES 


IF ( MODEL . EQ. 180 ) THEN 

DO 420 I - 1, NXI 
DO 421 J - 1, NYI 
DIFF ( I , J ) * -1000.0 
IF ( GA( I , J ) .NE. -1000.0 ) THEN 
DIFF( I , J) * GA ( I , J ) - GRA180 ( I , J) 
END IF 
42 1 CONTINUE 

420 CONTINUE 

ELSE IF ( M .EQ. 36 ) THEN 


DO 430 I - 1, NXI 
DO 431 J = 1, NYI 
DIFF ( I , J ) - -1000.0 
IF ( GA ( I , J ) .NE. -1000.0 ) THEN 
DIFF( I , J) - G A ( I , J ) - GRA036( I, J) 
END IF ' 

431 CONTINUE 

430 CONTINUE 

ELSE IF ( M .EQ. 300 ) THEN 


DO 440 1=1, NXI 

DO 441 J = 1, NYI 
DIFF ( I , J ) = -1000.0 

IF ( GA ( I , J ) .NE. -1000.0 ) THEN 
DIFF ( I , J ) « GA ( I , J ) - GRA3 00 ( I , J ) 
END IF 
441 CONTINUE 

440 CONTINUE 


WRITE ( OUP, 622 ) 

DO 450 I = Ml, M2 
WRITE ( OUP, 605 ) 
WRITE ( OUD , 606 ) 


( GRIDX ( I) , GRIDY (J) ,DIFF( I, J) , J=N1,N2 
( GRIDX ( I ) , GRIDY ( J ) , DIFF ( I , J ) , J=N1,N2 


) 

) 
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450 CONTINUE 


501 FORMAT ( IX , D 1 0 . 4 , 3D2 0 . 13) 

502 FORMAT( (4D20.13) ) 

503 FORMAT ( *6(2X,F10.3) ) ) 

504 FORMAT( ( 3 ( 2F7 . 2 , F12 . 4 ) ) ) 

505 FORMAT( ( 3X, 6511) ) 

600 FORMAT ( 1 H 1 ) 

601 TORMAT( 1H1,T 5, 'VALUES OF THE COVARIANCE TABLE USED', 

* ' IN ITERATION: ' , 13 , / 1H0, / ) 

602 FORMAT ( ( 1H , 5X, D 10 . 4 , 3 * 5X, D20 . 13 > ) ) 

603 FORMAT* 1H0,T2, 'GRID* ' ,12, ',' ,12, ')' ,2X, ' ELON =• ',F6.2,1X, 

* ' , GLAT =* ' ,F6 .2, 5X, 'GRAVITY ANOM - ',F10.4,' +- ',F6.2, 

* ', DELTA = ' , F6 . 2 , 5X, ' REC #: ',I3,2X,I8) 

604 FORMAT* 1H1,T5, 14, ' GRID POINTS HAVING NEGATIVE VARIANCE ',/, 1H0 ,/ , 

* T5, 'THERE ARE ',14,' POINTS GRIDDED OUT OF A NET ', 

* ' OF ',15,' GRID POINTS. ',/, 1H0, /,T8, 'THE RANGE or THE', 

* ' GRAVITY ANOMALY IN MGAL IS: ' , 2F12 . 4 , / , 1H0 ) 

605 FORMAT* ( 1H , 5* 2F7 . 2 , F12 . 4 ) ) ) 

606 FORMAT* ( 3 < 2F7 . 2 , F12 . 4 ) ) ) 

607 FORMAT* 1H1,T5, 'THE RANGE OF THE VARIANCE IN MGAL IS : ' , 2F12 . 4 , / 1H0 ) 

608 FORMAT* 1H1,T10, 'MODEL GRAVITY ANOMALY AND GEOID UNDULATION :',/ 1H0 ) 

609 FORMAT* ( 1H , T5 , 2F10 . 2 , 6F12 . 3 , 4X , 2F14 . 4 ) ) 

610 FORMAT* 1H1,T10, ' INPUT GEOID UNDULATION :',/, 1H0 ) 

611 FORMAT* 1H1,T10, ' INPUT GEOID STANDARD ERROR: ',/, 1H0 ) 

612 TORMAT* 1H1,T10, 'THE BLACK SEA MASK : ' , / , 1H0 ) 

613 FORMAT* * 1X,5(2F7.2,F12.4) ) ) 

614 FORMAT* 1H0, //,T8, 'THE RELATIVE RMS DIFFERENCES IS : ',F10.4,/, 

* T8, 'THE RMS RESPECT TO MODEL 300 IS : ',F10.4,/, 

* T8, 'THE RMS RESPECT TO MODEL 180 IS : ',F10.4,/, 

* T8 , 'THE RMS RESPECT TO MODEL 36 IS : ',F10.4) 

615 FORMAT* 1H1,T5, 14, ' GRID POINTS HAVING NEGATIVE VARIANCE ’,/, 1H0 ,/ , 

* T5, 'THERE ARE ',14,' POINTS GRIDDED OUT OF A NET ', 

* ' OF ',15,' GRID POINTS. ',/, 1H0, /,T8, 'THE RANGE OF THE', 

* ' GEOIDAL HEIGHT IN METER IS : ' , 2F12 . 4 , / , 1H0 ) 

616 FORMAT* 1H1,T5, 'THE RANGE OF THE VARIANCE IN CM IS : ' , 2F12 . 4 , / 1H0 ) 

617 FORMAT* 1H0,T2, 'GRID* ' ,12, ',' ,12, ')' ,2X, 'ELON » ',F6.2,lX, 

* ' , GLAT * ' ,F6.2,5X, 'GEOID HEIGHT » ',F10.4,' +- ' , F6 . 2 , 

* ', DELTA * ' ,F6.2,5X, 'REC #: ’,I3,2X,I8) 

618 FORMAT* 1H1 , ///,T30, 'OSU MODEL USED: ',13,////, 

T10, 'SUMMARY OF THE RMS Or THE ITERATIONS:',///, 


♦ 

T10 , 

1 RELATIVE 

’, T26, 

" RELATIVE 

* 

T40 , 

' 300 X 300 

T56, 

' 300 X 300 

★ 

T70, 

' 180 X 180 

', T86 , 

' 180 X 180 

♦ 

T100 , 

' 36 X 36 

' , Tl 16 , 

' 36 X 36 

* 

T2 , ' ITER' , 



♦ 

T10 , 

'GRAVITY ANOMALY’ 

, T26, 

'GEOID HEIGHT 


T40 , 'GRAVITY ANOMALY', T56, 'GEOID HEIGHT', 
T70, 'GRAVITY ANOMALY', T86, 'GEOID HEIGHT', 
T 100, 'GRAVITY ANOMALY ’, T1 16 ,' GEOID HEIGHT'/ 

( METER ) ' , 
* METER ) ' , 
( METER ) 

( METER ) ' / 


# 

T10 , ' 

(MGAL) 

T26, 

♦ 

T40 , ' 

(MGAL) 

T56 , 

* 

T70, ' 

(MGAL) 

', T86, 

★ 

T 1 0 0 , ' 

(MGAL) 

' ,T116, 

★ 

T2, ’ ' 



★ 

T10 , ' 


, T26, ' 

★ 

T4 0 , ' 

— 

, T56, ' 

★ 

T70, ' 


-' , T86 , ' 

* 

T100, ' 

— 

— ' , Tl 16 , ' 


621 FORMAT* 
* 

622 FORMAT* 
★ 

631 FORMAT* 
* 

632 FORMAT* 
* 

633 FORMAT* 

634 FORMAT* 


4 , T42 , F9 . 4,T57, F9 . 4 


( 1H0,/,T2,I3,T12,F9.4,T27,F9. 

T72,F9.4,T87,F9.4,T102,F9.4,T117,F9.4) ) 

1H 1 , T5 , ' THE DIFFERENCE BETWEEN THE PREDICTED AND ’, 

' THE MODEL GRAVITY ANOMALIES :',/, 1H0 ) 

1H ,T5, 'VALUES OF THE COVARIANCE TABLE USED’, 

' IN ITERATION: ’,13) 

1H ,T5, 'VALUES OF THE CORRESPONDING POWER SPECTRA USED' 
' IN ITERATION: '',13) 

IX, 110, 3D20. 13) 

1H 1 , T5 , 'VALUES OF THE CORRESPONDING POWER SPECTRA USED' 
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c 

c 


c 

c 

c 

c 


c 


c 

c 

c 


. * ' IN ITERATION: *,I3,/,1H0) 

635 FORMAT (( 1H , 5X, I 10 , 3 ( 5X, D20 . 13 ) ) ) 

999 FORMAT( 1H1,//,T10, '*♦** WRONG MODEL: ',15,' USED. 
* T20, 'APPROPRIATE MODELS ARE 36, 180 


STOP 

END 

SUBROUTINE POWSPC ( A, CV, X, 


P, LOT, LCTF, NFFT, 


.THIS SUBROUTINE DETERMINE THE POWER SPECTRUM Or A 
FUNCTION USING IMSL FFT SUBROUTINE 


REAL*8 
REAL* 8 
COMPLEX* 16 
INTEGER 


A ( LCTr ) , CV(LCT), P(LCTF), WK(NrFT) 
DI, DR, PMAX 
X ( LCTF ) 

IWK(NFFT) 


PMAX ■ -1.0D0 
ND2 » LCTF/ 2 + 1 


FIRST SYMMETRIZE THE COVARIANCE FUNCTION 


DO 100 I » 1, LCT-1 
A( I ) - CV( LCT+l-I ) 

A( LCT- 1+1 ) - CV(I) 

100 CONTINUE 
C 

c THEN SET THE END POINTS TO HAVE ZERO SLOPE 


****', //, 

AND 300 1 ' ) 


IWK, WK ) 
COVARIANCE 


A( 1) - A(2) 

A ( LCTF ) - A(2) 


CA LI* FFTRC ( A, LCTF, X, IWK, WK ) 
c DETERMINE THE REMAINING COEFFICIENTS 

c 

DO 110 1*2, ND2 

X( LCTF+2-I ) » DCONJG< X ( I ) ) 

110 CONTINUE 
C 

c DETERMINE THE POWER SPECTRUM 

C 


C 

C 

C 


C 


c 


DO 120 1*1, LCTF 

DR * DREAL ( X ( I) ) 

DI - DIMAG ( X( I ) ) 

P(I) * DSQRT ( DR*DR + 
120 CONTINUE 

..•.NORMALIZE THE COVARIANCE 

DO 130 I » 1, LCTF 

PMAX - DMAXl ( PMAX, P 
130 CONTINUE 

DO 140 1=1, LCTF 

P( I ) = P ( I ) /PMAX 
140 CONTINUE 


DI*DI ) 

FUNCTION BEFORE THE FOURIER 
I) ) 


TRANSFORM 


C 


RETURN 

END 

SUBROUTINE EMPCVF ( D,CTV, M, N , DXY , MN , G, DG, H , DH , CG,CH 

CHHNS , CHHEW , CHGNS , CHGEW , CGGNS , CGGEW ) 


REALM 
REALM 
REALM 
REAL *4 
REALM 
REALM 


CG<M,N), CH ( M, N ) 

G{M,N), DG ( M , N ) , H { M, N ) , DH ( M, N ) 
CHHNS ( 0 : MN- 1 ) , CHHEW ( 0 : MN- 1 ) 
CHGNS ( 0 : MN- 1 ) , CHGEW{ 0 :MN-1 ) 
CGGNS ( 0 : MN- 1 ) , CGGEW ( 0 : MN- 1 ) 

D ( MN ) , CTV { MN , 3 ) 


100 


non non 


REMOVE THE MODEL VALUE FROM THE DATA 
DO 100 I - 1, M 


►101 J - l, 

N 





IF ( DH(I,J) 

.NE. -1000.0 .AND. 

H( I, J) 

.NE. 

- 1000.0 ) 

THEN 

CH ( I , J ) - 

DH(I,J) - H ( I , J) 




END IF 

IF ( DG ( I , J ) 

.NE. -1000.0 .AND. 


.NE. 

- 1000.0 ) 

THEN 


CG(I,J) - DG(I,J) - G( I , J ) 
END IF 


101 CONTINUE 
100 CONTINUE 

PERFORM THE CONVOLUTION 

DO 200 K - 0, MN-1 
C 

SHH - 0. 

SHG - 0. 

SGG * 0. 

CHHNS(K) » 0. 

CHGNS(K) * 0. 

CGGNS(K) * 0. 

DO 210 J ■ 1, N 
DO 211 I ■ K+l, M 

IF ( CH< I-K, J) .NE. -1000. 0 .AND. CH( I , J) . NE . - 1 000 . 0 ) THEN 
SHH - SHH + CH< I-K, J)*CH( I, J) 

END IF 

IF ( CH(I-K, J) .NE. -1000.0 .AND. CG( I, J) .NE. -1000 . 0 ) THEN 
SHG * SHG + CH( I-K, J) *CG( I , J) 

END IF 

IF ( CG(I-K, J) .NE. -1000.0 .AND. CG( I , J) . NE . -1000 . 0 ) THEN 
SGG * SGG + CG ( I-K, J ) *CG ( I , J ) 

END IF 
2 1 1 CONTINUE 

210 CONTINUE 

Cl - ( (M-K) *N) 

CHHNS(K) - SHH / Cl 
CHGNS(K) - SHG / Cl 

CGGNS(K) =» SGG / Cl 

C 

SHH = 0. 

SHG =■ 0. 

SGG - 0. 

CHHEW(K) - 0. 

CHGEW(K) » 0. 

CGGEW(K) =* 0. 

DO 220 I = 1, M 
DO 221 J - K+l, N 

IF ( CH( I, J-K) .NE.-1000. 0 .AND. CH ( I , J) . NE . - 1000 . 0 ) THEN 
SHH =* SHH + CH( I, J-K) *CH( I, J) 

END IF 

IF ( CH< I, J-K) .NE. -1000.0 .AND. CG( I , J) . NE . - 1 000 . 0 ) THEN 
SHG = SHG + CH ( I , J-K ) *CG ( I , J ) 

END IF 

IF ( CG( I, J-K) .NE. -1000.0 .AND. CG< I , J) . NE . - 1 000 . 0 ) THEN 
SGG » SGG + CG( I, J-K) *CG( I, J) 

END IF 
221 CONTINUE 

220 CONTINUE 

C2 = ( M* ( N-K ) ) 

CHHEW(K) = SHH / C2 

CHGEW(K) = SHG i C2 

CGGEW(K) = SGG / C2 

C 

D(K+1) = DBLE ( K * DXY ) 

CTV ( K+ 1,1) = DBLE (< C1*CHHNS(K) + C2*CHHEW(K) ) / ( Cl + C2 )) 

CTV ( K+ 1,2) = DBLE (( C1*CHGNS(K) + C2*CHGEW(K) ) / ( Cl + C2 ) ) 

CTV ( K+ 1,3) = DBLE (( Cl*CGGNS(K) + C2*CGGEW(K) ) / ( Cl + C2 )) 
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c 


c 


c 


c 

c 

c 


200 CONTINUE 


RETURN 

END 

SUBROUTINE PREDIC ( IT, PVAL,VAR, EX , GY , GH , SGY,CGY, IREC, 

* SE, X1,Y1,SY1,CYI, 

* CM, CMS , CMI , JREC , IDIAG, CP, SRC, AP, 

* SD , CTV , CVUU , SCUU , CVUG , SCUG , CVGG , SCGG , 

* LCT, DTR, IER, NIER ) 


REAL *8 
REAL *8 
REAL *8 
REAL *8 
REAL *8 
REAL *8 
REAL*8 
REAL *8 
REAL* 8 
REAL *8 
REAL *8 
INTEGER 


SE( IREC) 

EX(IREC), GY ( IREC ) , GH(IREC), SGY(IREC), 
CP (IREC), SRC (IREC), AP(IREC) 

CM ( JREC ) , CMI ( JREC ) , CMS(JREC) 

SD(LCT) , CTV ( LCT , 3 ) 

CVUU ( LCT ) , SCUU ( LCT- 1,3) 

CVUG (LCT), SCUG (LCT- 1,3) 

CVGG(LCT), SCGG < LCT- 1,3) 

Dl, D2 , CSD, DSD, VDSD, DTR 
XI, Yl, SY1, CY1, PVAL, VAR 
Rl, R2 
IDIAG ( IREC) 


CGY ( IREC) 


RESET THE ERROR COUNTER 


IER - 0 
C 

IF ( IT. EQ. 1 ) THEN 
Rl - CTV ( 1,1) 

R2 - CTV ( 1,3) 

ELSE IF ( IT. EQ. 2 ) THEN 
Rl - CTV( 1, 3) 

R2 = CTV ( 1 , 1) 

END IF 
C 

c CLEAR THE COVARIANCE MATRIX ARRAY 


C 


C 

c 

c 


c 

c 

c 


c 


CALL DCLR( CM, JREC ) 
CALL DCLR( CMI, JREC ) 
CALL DCLR< CMS, JREC ) 


DO 100 1*1, 

CM( IDIAG ( I ) ) 
100 CONTINUE 


IREC 

* Rl + SE( I) *SE( I ) 


....EVALUATE ALL THE SINE AND COSINE TERMS 

DO 200 1*1, IREC 

SGY(I) - DSIN ( GY ( I ) ) 

CGY ( I ) * DCOS ( GY ( I ) ) 

200 CONTINUE 


CONSTRUCT THE OFF-DIAGONAL TERMS OF THE SIGNAL COVARIANCE 

DO 3 00 J * 2, IREC 
DO 301 1=1, J-l 


MATRIX 


C 


CSD = SGY ( I ) *SGY ( J ) + CGY ( I ) 
DSD » DACOS(CSD) /DTR 
IF ( IT. EQ. 1 ) THEN 

CALL ICSEVU ( SD, CVUU, LCT, 
ELSE IF ( IT. EQ. 2 ) THEN 

CALL ICSEVU ( SD, CVGG, LCT, 
END IF 

K = < J*(J-1) )/2 + I 
CM( K ) = VDSD 


♦CGY ( J ) *DCOS ( EX ( I ) -EX ( J 


SCUU, LCT- 1, 
SCGG, LCT- 1, 


DSD, VDSD, 
DSD, VDSD, 


) 

1 , 

1 , 


301 CONTINUE 
300 CONTINUE 


IER ) 
IER ) 
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c 

C DUPLICATE THE COVARIANCE MATRIX FOR STORAGE 

C 

DO 310 1*1, JREC 

CMS ( I ) » CM( I ) 

310 CONTINUE 


C 

C INVERT THE COVARIANCE MATRIX 

C 

CALL LINV1P ( CMS , I REC , CMI , IDGT, Dl,D2, IER ) 

C CALL SOLVE ( CMS, SRC, IDIAG, IREC, .TRUE,, .FALSE. ) 

C CALL SYMINV ( CMS, SRC, IDIAG, IREC, CMI ) 

C 

IF ( IER. EQ. 0 ) THEN 

C 

C CONSTRUCT THE COVARIANCE VECTOR 


C 

DO 330 I - 1, IREC 
C 

CSD * SY 1*SGY ( I ) + CYl*CGY(I)*DCOS(Xl-EX(I) ) 

DSD * DACOS( CSD) /DTR 

CALL ICSEVU ( SD , CVUG , LCT , SCUG,LCT-1, DSD, VDSD, 
CP(I) - VDSD 
C 


330 CONTINUE 
C 

C DETERMINE THE PREDICTED VALUE AT A GRID POINT 

C 


CALL VMULSF < CMI, IREC, GH,1,IREC, SRC, IREC ) 

C CALL DCLR ( SRC, IREC ) 

C CALL VSMXVT ( CMI, GH, SRC, IDIAG, IREC ) 

PVAL * DOT ( CP, SRC, IREC ) 

C 

C DETERMINE THE VARIANCE OF THE PREDICTED VALUE 

C 

CALL VMULSF ( CMI, IREC, CP, 1, IREC, AP,IREC ) 

CALL VMULSF ( CM, IREC, AP,1,IREC, SRC, IREC ) 

C CALL DCLR ( AP, IREC ) 

C CALL DCLR ( SRC, IREC ) 

C CALL VSMXVT ( CMI, CP, AP, IDIAG, IREC ) 

C CALL VSMXVT ( CM, AP, SRC, IDIAG, IREC ) 

VAR - R2 - 2 . 0D0*DOT ( AP , CP , IREC ) + DOT ( AP , SRC , IREC ) 

C 

IF ( VAR. GT. 0.0 ) THEN 
VAR * DSQRT(VAR) 

ELSE 

IER * 1 

NIER = NIER + 1 
END IF 
C 

END IF 
C 

RETURN 

END 

FUNCTION DOT ( A , B , N ) 

C 

C Compute the dot product of the two N-vectors A and B. 

C 


C 

C 


C 


INTEGER N, I 
REAL* 8 DOT 

REAL* 8 A ( 1 ) , B ( 1 ) 

DOT=0 . 0 


DO 100 1=1, N 

100 DOT=DOT + A(I)*B(I) 

RETURN 

END 

SUBROUTINE SYMINV ( 


A, 


B, 


IDIAG, NEQ, AINV ) 


, IER ) 
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IMPLICIT REAL *8 (A-H,0-Z) 

DIMENSION A(l), B<1), AINV(l), IDIAG(l) 

K * 0 

DO 100 1-1, NEQ 

DO 101 L * 1 , NEQ 
B(L) - 0.0D0 
101 CONTINUE 

B(I) - 1.0D0 

CALL SOLVE ( A, B, ID I AG, NEQ, .FALSE., .TRUE. ) 

DO 110 J - 1, I 
K * K + 1 
AINV(K) - B(J) 

110 CONTINUE 
100 CONTINUE 

RETURN 

END 

SUBROUTINE SOLVE ( A, B, IDIAG, NEQ, FACT, BACK) 

Compute the U**T * D * U factorization of the symmetric matrix 
stored m A, if FACT » TRUE; and solve A * X ■ B if BACK * true. 


A 


B 


IDIAG 

NEQ 

FACT 

BACK 


Contains the compacted-column form of the upper triangular 
part of the coefficient matrix. After factorization, it 
contains D and U. 

Right-hand-side vector. After backsubstitution , it 
contains the solution. 

Addresses of the diagonal terms in A. 

Number of equations 

If FACT * TRUE, then factor A; otherwise do not factor A. 
If BACK * TRUE, reduce B and backsubstitute; otherwise 
do not solve the equations. 


IMPLICIT REAL *8 (A-H,0-Z) 

LOGICAL FACT, BACK 
DIMENSION A( 1 ) , 8 ( 1 ) , IDIAG ( 1 ) 


Factor A, reduce B 
JR * 0 

DO 400 J - 1, NEQ 
JD - IDIAG ( J ) 

JH - JD - JR 
IS * J - JH + 2 
C 

IF (JH .LT. 2) GOTO 390 
C 

IF (FACT) THEN 
C 

IF (JH .GT. 2) THEN 


C. . 

Reduce column J rows is 
K * JR + 2 
ID « IDIAG ( IS - 1) 

to J-l: 

c 



DO 100 I =» IS, J-l 



IR =* ID 



ID = IDIAG ( I ) 

IH = MIN ( ID-IR-1 , I- 
IF (IH.GT.0) A(K) = 
K = K + 1 

IS+1) 
A(K) - 

100 

CONTINUE 

ENDIF 



do not divide by row diagonal 


DOT ( A ( K- I H ) ,A( ID-IH) , IH) 
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Divide by row diagonal, and reduce diagonal term in column J 

IR » JR + 1 
K » J - JD 
DO 200 I » IR, JD-1 
ID ■ IDIAG ( K+I ) 

IF (A(ID) .EQ.0.0) GOTO 200 
D - -A( I ) 

A( I ) - A( I ) /A( ID ) 

A(JD) - A( JD ) + D*A( I ) 

200 CONTINUE 

ENDIF 

. Reduce RHS 

IF (BACK) B(J) » B ( J ) - DOT ( A( JR+1 ) , B( IS— 1 ) , JH-1 ) 

390 JR - JD 
400 CONTINUE 


IF (.NOT. BACK) RETURN 


Divide by diagonal pivots 

DO 700 I m 1 , NEQ 
ID - IDIAG ( I ) 

IF (A( ID) .NE. 0.0) B ( I ) - B(I)/A(ID) 
700 CONTINUE 

Backs ubstitute 

J - NEQ 

JD - IDIAG ( J) 

801 D » -B(J) 

J - J - 1 

IF (J.LE.0) RETURN 

JR - IDIAG ( J) 

IF ( JD-JR.GT. 1 ) THEN 
IS * J - JD + JR + 2 
K - JR - IS + 1 
DO 810 I - IS, J 
810 B ( I ) * B ( I ) + A(I+K)*D 

ENDIF 

JD * JR 
GOTO 801 


END 

SUBROUTINE VSMXVT ( A, B, C, JDIAG, NEQ ) 

THIS SUBROUTINE FORMS C * C + A*B WHERE A IS A SYMMETRIC MATRIX 
STORED IN PROFILE FORM, B, C ARE VECTORS, AND JDIAG LOCATES THE 
DIAGONALS IN A. 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A(l), B(l), C(l), JDIAG(l) 

JS = 1 

DO 200 J = 1, NEQ 
JD = JDIAG ( J ) 

IF ( JS. LE. JD ) THEN 
BJ = B ( J ) 

AB = A( JD) *BJ 
IF ( JS. NE. JD ) THEN 
JB = J - JD 
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c 


100 


200 


JE - JD - 1 

DO 100 JJ ■ JS, JE 

AB * AB + A( JJ) *B( JJ+JB) 

C ( JJ+JB ) m C( JJ+JB) + A( JJ) *BJ 
CONTINUE 
END IF 

C(J) - C(J) + AB 
END IF 
JS - JD + 1 
CONTINUE 


RETURN 

END 

SUBROUTINE DCLR ( A , NA ) 
C 

INTEGER NA, I 
REAL* 8 A( 1 ) 

C 

DO 100 I - 1, NA 
A(I) * 0.0D0 
100 CONTINUE 
C 


BLACK SEA GRIDDED COLLOCATION GEOIDAL HEIGHT DATA 


RETURN 

END 

//* 

// EXEC LINKGOV, 

//* 

//SYSLIB DD DSN 
//* DD DSN 

//* 

//* 

//* 

//* 

//* 

//* * + 

//* 

//Go.rrosrooi oo 
//Go.rT06rooi dd 
//Go.rToerooi dd 
//GO.rT09F001 DD 
//GO.FT10F001 DD 
/ /GO. FT11F001 DD 
/ /GO. FT12F00 1 DD 
/ /GO. FT13F001 DD 
//GO.FT14F001 DD 
//GO. FT15F00 1 DD 
/ /GO. FT16F00 1 DD 
/ /GO. FT17F001 DD 
/ /GO.FT10FOO1 DD 
//GO.FT19F001 DD 
//GO.FT20F001 DD 
//GO.FT21F001 DD 
/ /GO. FT22F001 DD 
//* 

// EXEC NOTIFYTS 


REGION-5000K 

SYS2 . IMSLD, DISP-SHR 
-SYS2 . IMSLS, DISP-SHR 


DSN-ZMAYA 

SYSOOT-* 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA 

DSN-ZMAYA. 

DSN-ZMAYA. 

DSN-ZMAYA 


ALTIM. DATA < TABLE 1 ) , DISP-SHR 

. BLACKC4H . GRID, DISP-SHR 
. BLACKC4H . COVM , DISP-SHR 
. ALTIM.DATA(BRH1 80 ), DISP-SHR 
. ALTIM. DATA( BRG180) , DISP-SHR 
. ALTIM . DATA ( BRH 36) , DISP-SHR 
.ALTIM. DATA (BRG3 6) , DISP-SHR 
.ALTIM. DATA (BRH3 00 ), DISP-SHR 
• ALTIM. DATA( BRG300 j , DISP-SHR 
. BLACK . MASK , DISP-SHR 
•ALTIM. DATA ( BSGAB0R) , DISP-SHR 
ALTIM. COVM ( BSGAB0R ) , DISP-SHR 
ALTIM. DATA ( BSGUB0R) , DISP-SHR 
ALTIM. COVM (BSGUB OR ), DISP-SHR 
ALTIM. DIFF ( BSGUB0R) , DISP-SHR 
ALTIM. CVTB( BOR) , DISP-SHR 


* 


♦ 
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A listing Of PROGRAM JORDAN 


//ZMAYABJC JOB ( GO 109 , 360 , 2 ) , AYAU, TIME>( 00 , 30 ) , CLASS-O, MSGCLASS 
// EXEC FORTVC 
//SYSIN DD * 

C 

PARAMETER ( SN - 1.27 ) 

PARAMETER ( SG - 11.20 ) 

PARAMETER ( D » 0.38 ) 

PARAMETER ( DR - 0.25 ) 

PARAMETER ( RMIN « 0 . 0 ) 

PARAMETER ( RMAX =» 10.0 ) 

PARAMETER ( SN2 - SN*SN ) 

PARAMETER ( SG2 = SG*SG ) 

PARAMETER ( L - RMAX/DR + 1 ) 


C 

C 


c 


REAL* 8 
REAL *8 
REAL* 8 
REAL *8 
INTEGER 


R(L), CNN ( L) , CNG(L), CGG ( L ) 

RD, RDS, RD2 , EXPRD, BIO, BI1, BKO, BK1 
DSN2, DSG2 , DD, DDR, CNST 
MMBSI0, MMBSI1, MMBSKO, MMBSK1 
OUP, OUC 


DATA 


OUP / 6 /, OUC / 8 / 


DSN2 - DBLE( SN2 ) 

DSG2 ■ DBLE ( SG2 ) 

DDR » DBLE (DR) 

DD - DBLE(D) 

CNST - DBLE(2.*SN*SG)/DSQRT(«.0D0) 


R( 1) - DBLE (RMIN) 

DO 100 I » 2, L 

R(I) - R(I-l) + DDR 
100 CONTINUE 


CNN ( 1) * DSN2 
CGG(l) - DSG2 
CNG( 1 ) - CNST 
DO 200 I - 2, L 


RD - R(I)/DD 
RDS * RD*RD 
RD2 = RD/2.0D0 
EXPRD = DEXP(-RD) 

BIO * MMBSI0 ( 1 , RD2 , IER) 
BI1 - MMBSI 1(1, RD2 , IER ) 
BKO - MMBSKO ( 1,RD2, IER) 
BK1 =■ MMBSK 1(1, RD2 , IER ) 


X 


CNN ( I ) - DSN2 * ( 
C 

CGG(I) - DSG2 * ( 

c 

CNG(I) = CNST * ( 
♦ + 
C 

200 CONTINUE 
C 


1.0D0 + R D + RDS/3.0D0 ) * EXPRD 

1.0D0 + RD - RDS/2.0D0 ) * EXPRD 

RD2* { 1 . 0D0- ( RDS/ 2 . 0D0 ) )*(BIO*BK1-BI1*BKO) 
( RDS/4 . 0D0 ) * ( BI0*BK0+BI 1*BKI ) ) 


C 


c 


WRITE ( OUP, 600 ) 

WRITE ( OUP, 601 ) 

WRITE { OUC, 801 ) 


( R( I ) , CNN { I ) , CNG ( I ) , CGG ( I ) , 1 = 1 ,L) 
(R{ I) , CNN ( I ) , CNG ( I ) , CGG ( I ) , I=i,L) 


6 00 ^FORMAT (1H1,//,T5,' THE COVARIANCE FUNCTION C(NN), C ( NG ) 

601 FORMAT ( ( 1H , 5X, Dl 0 . 4 , 3 ( 5X , D2 0 . 13 ) ) ) 

801 FORMAT ( { 1X,D10.4,3D2 0. 13) ) 


& C { GG ) 


STOP 

END 
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//* 

// EXEC LINKGOV, REGION*5000K 

//* 

//SYSLIB DD DSN-SYS2.IMSLD,DISP-SHR 
//* DD DSN-SYS2.IMSLS,DISP-SHR 

//GO. FT06F001 DD SYSOUT** 

//GO.FT08F001 DD DSN-2MAYA. ALTIM. DATA( BJ1 80 ) , DISP-SHR 
// EXEC NOTIFYTS 
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